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INTRODUCTION 


A  research  program  of  basic  experimental  and  computational  studies  on  the  unsteady 
flow  generated  by  a  pitching  airfoil  was  initiated  in  October  1986  under  grant  no.  AFOSR- 
86-0243.  This  report  is  being  written  to  describe  the  work  carried  out  during  the  period 
October  1986  -  October  1988.  It  is  of  interest  to  note  that  the  work  is  being  continued 
under  a  separate  contract  AFOSR  F49620-89-C-0014. 

The  present  study  has  the  following  objectives  to  the  general  problem  of  unsteady 
flow  past  a  NACA  0012  airfoil  undergoing  large  amplitude  incidence  variation. 

1.  Investigate  the  unsteady  flow  structure  generated  by  an  impulsively  started  airfoil 
from  rest,  at  different  angles  of  attack. 

2.  Investigate  the  transient  flow  characteristics  of  an  NACA  0012  airfoil  undergoing  a 
stepwise  varying  angle  of  attack. 

3.  Investigate  the  unsteady  flow  structure  of  an  accelerating  NACA  0012  airfoil  at  dif¬ 
ferent  fixed  angles  of  attack. 

4.  Investigate  the  transient  flow  characteristics  of  an  accelerating  airfoil  undergoing  a 
stepwise  incidence  variation. 

In  all  these  studies  attention  was  directed  to  the  basic  understanding  of  the  unsteady 
flow  phenomena.  A  unique  experimental  technique,  known  as  Particle  Image  Displacement 
Velocimetry,  developed  in  our  laboratory,  is  successfully  implemented  to  study,  in  detail, 
the  unsteady  large  scale  vortical  motions  that  occur  in  these  flows.  A  parallel  effort  was 
devoted  to  study  the  experimentally  observed  features  using  numerical  simulations.  In  an 
attempt  to  expand  the  parameter  range,  in  consultation  with  AFOSR,  a  new  larger  towing 
tank  facility  has  been  designed  and  constructed. 

In  the  following,  the  status  of  the  research  effort  is  given  along  with  some  of  the  most 
important  conclusions  arrived  from  this  study.  The  details  of  various  investigations  are 
given  in  reports  and  papers,  which  are  included  here  as  appendices. 


STATUS  OF  THE  RESEARCH  EFFORT 


In  the  following,  major  conclusions  on  the  various  facets  of  the  research  program  are 
given. 

a)  Impulsively  stated  airfoil. 

A  systematic  investigation  has  been  carried  out  to  understand  the  basic  flow  structure 
generated  by  an  impulsively  started  NACA  0012  airfoil,  of  finite  aspect  ratio,  at  different 
angles  of  attack,  and  at  a  fixed  Reynolds  number  of  1400. 

A  novel  experimental  technique  known  as  “Particle  Image  Displacement  Velocimetry” 
was  used  to  measure  the  instantaneous  two  dimensional  velocity  field.  The  velocity  field 
was  measured  with  sufficient  accuracy  cind  spatial  resolution  that  the  vorticity  field  and 
pressure  field  can  be  computed  accurately,  a  unique  capability  of  the  technique.  The 
detailed  description  of  this  technique  is  given  in  the  Appendix  I. 

A  parallel  computational  study  was  c-jnducted  to  augment  the  above  mentioned  ex¬ 
perimental  study.  In  this  study,  a  new  algorithm  was  developed  to  solve  the  Navier-Stokes 
equations  using  the  discrete  vortex  method.  The  new  fast  velocity  summation  algorithm 
enables  the  flow  to  be  computed  with  much  more  resolution  than  previously  possible  in 
vortex  methods.  The  details  of  the  algorithm  is  given  in  Appendix  II. 

The  main  features  of  the  unsteady  large  scale  separated  flow  about  an  impulsively 
started  airfoil  are  as  follows; 

The  multiple  exposure  photographs  of  the  flow  field  about  the  airfoil  at  lO"  or  less 
incidence  showed  that  the  flow  is  well  behaved  and  attached  to  the  airfoil  over  the  period 
of  observation.  However,  at  large  angles  of  attack  a  >  20",  the  flow  separates  on  the  upper 
surface  of  the  airfoil  and  generates  large  scale  vortices.  The  following  scenario  develops 
in  time  on  the  upper  surface  of  the  airfoil.  At  the  start  of  the  airfoil,  a  starting  vortex 
develops  at  the  trailing  edge  and  is  carried  away  downstream  of  the  airfoil.  Concomitant 
with  is  the  generation  of  a  separation  bubble  at  the  leading  edge  of  the  airfoil.  At  a  later 
time,  the  separation  bubble  grows  into  an  i,solated  primary  vortex  with  secondary  vortices 
following  behind  it.  This  multiple  vortex  structure  continues  to  grow  together  and  move 
along  the  upper  surface  until  it  reaches  the  trailing  edge.  .\t  this  point,  the  primary  vortex 
in  duces  a  vortex  at  the  trailing  edge  with  the  sense  opposite  to  that  of  the  primary  vortex 
Finally,  the  primary  vortex  and  the  induced  trailing  edge  vortex  interact  and  generate  a 
complex  flow  field.  However,  for  finite  aspect  ratio  airfoils  or  wings,  a  different  type  of  flow 
field  seems  to  emerge  at  later  stages  of  development.  The  various  events  described  above 
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occur  at  different  times,  depending  upon  the  zmgle  of  attack  and  free  stream  Reynolds 
nuniher. 

Typical  PIDV  measurements  of  the  instantaneous  velocity  field,  at  different  times,  for 
an  airfoil  at  a  =  30“  are  shown  in  figure  l.The  aspect  ratio  of  the  wing  was  about  3.  The 
airfoil  travels  from  right  to  left.  The  data  is  presented  in  the  body  fixed  reference  frame. 
The  length  of  the  velocity  vector  corresponds  to  its  magnitude.  The  dimensionless  time  t* 
is  defined  as  where  U  is  the  free  stream  velocity  2uid  C  is  the  airfoil  chord.  The  starting 
vortex  (at  the  far  right  of  the  picture)  and  the  intial  separation  bubble  at  the  leading  edge 
can  be  seen  in  the  figure  corresponding  to  t*  =  0.68.  The  primary  vortex  with  secondary 
vortices  following  behind  it  can  be  seen  in  the  figure  at  t*  =  2.02.  The  trailing  edge  vortex 
can  be  seen  in  the  figure  at  t*  =  3.02  At  t*  >  3,  the  primary  vortex  abruptly  moves  away 
from  the  upper  surface  leaving  behind  a  vortex  sheet  type  like  structure.  Such  a  behavior 
is  attributed  to  the  interference  of  tip  vortices  which  eire  generated  due  to  the  finite  aspect 
ratio  of  the  airfoil.  At  later  times,  for  example  at  i*  ~  4.85,  the  tip  vortices  interact  with 
the  separated  flow  on  the  upper  surface  and  generate  a  complicated  three-dimensional  flow 
field.  The  nature  of  this  interaction  and  parameters  that  govern  such  a  flow  field  is  not 
yet  known.  The  current  experiments  with  a  larger  aspect  ratio  (10)  airfoil  will  help  us  in 
the  interpretation  of  these  results. 

Typical  two-dimensional  computational  results  from  random-walk  vortex  simulations 
of  the  full  Navier-Stokes  equations  are  shown  in  figure  2.  The  angle  of  attack  and  the 
Reynolds  number  axe  the  same  as  those  in  the  experiment;  the  results  of  which  are  shown 
in  figure  1.  The  stream  line  pattern,  along  with  vorticity,  which  is  represented  in  bit¬ 
mapped  graphics  as  half  tones  are  shown  in  the  figure.  Elxcept  for  the  effect  of  the  finite 
aspect  ratio  of  the  airfoil,  the  stream  line  pattern  looks  very  much  similar  to  those  found 
in  figure  1.  To  further  evaluate  these  results,  the  locus  of  the  primary  vortex  as  it  develops 
in  time  is  shown  in  figure  3.  The  computational  results  agree  well  with  the  experiment  for 
t*  <  2.  Beyond  t*  =  2,  it  is  expected  that  the  experimental  flow  field  was  inflnenred  by  the 
tip  vortices  making  it  to  be  three-dimensional.  The  coefficients  of  lift  and  drag  as  obtained 
from  the  computations  are  shown  in  figure  4.  As  expected,  the  coefficient  of  lift  increases 
with  t"  up  to  a  point  where  the  primary  vortex  is  attached  to  the  upper  surface.  For 
later  times,  where  the  primary  vortex  leaves  the  upper  surface,  the  coefficient  of  lift  drops 
significantly.  In  order  to  have  a  proper  comparisons,  the  experiment  is  being  conducted 
at  this  time  with  a  larger  aspect  ratio  airfoil,  where  the  flow  can  remain  two-dimensiona' 
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beyond  t*  >  2. 

Similar  observaiions  are  also  madi  for  Mte  airfoil  at.  o  —  20",  atid  llio  results  are 
shown  in  figures  5  and  6.  Figure  5  shows  typical  flow  field  measurements  at  four  different 
times.  The  drop  out  regions  are  a  result  of  the  blockage  of  the  laser  light  by  the  airfoil. 
The  new  optical  set  up  explained  later  will  avoid  such  regions.  The  various  stages  of  the 
development  of  the  separated  flow,  as  explained  above,  can  be  seen  in  these  instantaneous 
flow  fields.  The  dynamics  of  these  vortices  and  their  role  on  the  surface  pressures  are  being 
investigated  at  this  time.  The  corresponding  computational  results  are  shown  in  figure  6. 
From  these  results,  it  appears  that  the  numerical  simulation  depicts  clearly  the  various 
stages  of  the  flow  development  observed  in  the  experiment. 

b)  Pitching  Airfoil. 

The  experimental  facility  to  investigate  the  transient  flow  characteristics  of  an  airfoil 
with  stepwise  varying  angle  of  attack  has  been  built  and  the  details  are  given  in  the  next 
section.  However,  numerical  simulations  have  been  carried  out  on  a  typical  case  of  a  NACA 
0012  airfoil  undergoing  a  rapid  variation  in  its  angle  of  attack  as  shown  in  figure  7,  In 
the  simulation,  the  airfoil  initially  translates  at  zero  angle  of  attack  during  a  time  interval 
0  <  t*  <  1;  in  the  interval  1  <  t*  <  2.3  the  angle  of  attack  increases  linearly  from  0“  to 
30";  beyond  this  time  the  angle  of  attack  is  held  fixed  at  30°.  The  numerical  procedure  is 
summarized  in  appendix  III.  The  results  of  this  numerical  simulation  are  shown  in  figure 
8.  The  results  clearly  shows  the  various  stages  of  the  transient  flow  development  and 
associated  vortical  flow  field.  These  results  are  still  being  analyzed. 

c)  Towing  Tank  Facility. 

A  computer  controlled  towing  tank  facility  has  been  designed,  and  constructed.  A 
schematic  of  the  facility  is  shown  in  figure  9.  The  various  components  of  the  facility  are 
indicated  in  the  figure.  The  significant  improvement  in  this  facility,  over  other  traditional 
towing  tanks,  is  the  ability  to  vary  the  velocity  rapidly  without  encountering  any  vibration 
of  the  model.  Which  is  accomplished  by  the  use  of  "Anorail”  linear  motor  system  controlled 
by  an  intelligent  axis  controller  with  a  speed  range  from  3mm/sec  to  0.3m/sec.  Impulsive 
or  continuous  airfoil  pitching  motion  is  controlled  by  ’’Klinger”  stepping  motor  system.  A 
new  optical  arrangement  has  been  incorporated  in  to  this  facility,  which  has  the  capability 
of  steering  the  laser  sheet  in  several  directions.  This  capability  will  enable  us  to  illuminate 
the  entire  flow  field  around  the  airfoil. 

The  laser  source  and  the  camera  are  synchronously  controlled  by  an  electronic  system 
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which  is  activated  by  a  computer.  A  schematic  of  the  arrangement  is  shown  in  figure  10. 
The  entire  operation  i<?  monitored  by  rtF.C  Vax  Station  TT  rompnter. 


5 


Figure  2.  Two-dimensional  computational  results  of  the  flow  field 
about  an  NACA  0012  airfoil.  Re  =  1400. 
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variation  of  the  primary  vortex  location  with  time 
3  of  attack  =30  degrees.  Re  =  1400. 


Figure  6.  Numerical  simulation  of  the  fiow  past  an  impulsively  started  airfoil  at  =20°j 
Re=  lAOO;  a)  t*=0.5;  b)t*=1.0;  c)t*=2.0:  d)t*=3.0;  e)t*=A.O;  f)t*=5.0. 
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Figure  7. 


The  variation  of  the  angle  of  attack  with  time  for 
the  computational  study. 
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(b)  Front  view 


schematic  of  the  towing  tank  arrangement 
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Abstract 

A  novel  experimental  technique, 
commonly  referred  to  as  Laser  Speckle 
Velocimetry  (LSV)  or  Particle  Image 
Displacement  Velocimetry  (PIDV),  is  devel¬ 
oped  for  the  measurement  of  instantaneous 
velocity  fields  in  unsteady  and  steady 
flows.  The  main  advantage  of  this  tech¬ 
nique  is  that  the  velocity  field  is 
measured  with  sufficient  accuracy  so  that 
the  distribution  of  vorticity  can  be  cal¬ 
culated  with  accuracy. 

The  PIDV  technique,  whicli  is  ideally 
suited  for  the  study  of  unsteady  separated 
flows,  has  been  utilized  to  measure  the 
development  of  the  separated  flow  field 
generated  by  an  high  angle  of  attack 
(o=30°)  NACA  0012  airfoil  started  impul¬ 
sively  from  rest. 

1.  Introduction 

For  the  solution  of  many  problems  that 
occur  in  high  angle  of  attack  aero¬ 
dynamics,  it  is  a  necessary  to  have  a 
thorough  understanding  of  the  behavior  of 
unsteady  separated  flows.  Although  much 
progress  has  been  made  in  predicting  the 
steady  flow  phenomenon  with  the  use  of 
numerical  methods,  it  is  still  difficult 
to  predict  unsteady  flows  which  contain 
flow  separation.  The  difficulty  mainly 
arises  from  the  fact  that  these  flows  are 
extremely  complex  and  are  not  amenable  to 
standard  experimental  and  numerical 
techniques.  In  view  of  this,  a  novel 
experimental  technique  is  being  developed 
for  the  measurement  of  instantaneous 
velocity  fields  in  unsteady  and  steady 
fluid  flows.  This  paper  provides  a 
description  of  this  technique  along  with 
its  successful  application  to  the  study  of 
an  unsteady  separated  flow  generated  by  ar 
high  angle  of  attack  airfoil. 

In  an  unsteady  flow,  a  single  photo¬ 
graph  of  the  flow  pattern  at  a  given 
instant  does  not  generally  provide  any 
meaningful  information.  In  order  to 
understand  the  unsteady  flow  phenomena,  it 
is  necessary  to  obtain  both  spatial  and 
temporal  information  of  the  entire  flow 
field.  With  this  in  mind,  optical 
techniques  have  been  widely  used  to 
observe  and  measure  properties  of  flow 
fields  such  as  velocities  and  densities. 
Many  of  these  techniques  are  qualitative 
in  nature,  but  of  great  value  in  guiding 
intuition  and  in  suggesting  ways  to 
investigate  the  problem  by  quantitative 
means.  An  admirable  review  of  the  modem 
optical  techniques  in  fluid  mechanics  is 
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given  by  Laucerborn  and  Vogel',  and  the 
reader  is  referred  to  this  article  for 
details.  A  ipiant i tat i ve  flow  visualiza¬ 
tion  technique  would  be  very  helpful  in 
the  study  of  these  flow  fields.  Attempts 
to  accom,  ..sh  this  task  by  tracing  the 
streaklines  of  injected  particles'' '  ' 
usually  can  not  provide  the  spatial 
resolution  that  is  required  and  a  large 
amount  of  labor  is  necessary  to  reduce  the 
data. 

In  unsteady  separated  flows,  as  will 
be  shown  later,  it  is  often  desirable  to 
obtain  the  vorticity  field,  in  addition  to 
tlie  velocity  field.  However,  measurement 
of  the  vorticity  exceed  the  present  exper¬ 
imental  capability.  This  difficulty 
arises  from  the  fact  that  vorticity  is  a 
quantity  defined  in  terms  of  local  veloc¬ 
ity  gradients.  In  contrast,  the  currently 
available  measurement  techniques,  such  as 
hot-wire  anemometry  or  laser  velocimetry, 
are  sensitive  only  to  the  local  velocity. 
Hence,  measurements  must  be  made  over 
several  points  and  the  resulting  velocity 
components  are  then  analyzed  by  finite 
difference  schemes.  However,  the  errors 
produced  by  the  necessary  differentiations 
limit  the  accuracy  and  spectral  range.  In 
add ition,  the  spatial  resolution  of  this 
metliod  is  oftt.n  not  sufficient  to  measure 
small-scale  tluid  motions  or  rapidly 
changing  velocity  gradients.  As  a  conse¬ 
quence,  the  measured  vorticity  field  is  a 
type  of  spatially  averaged  estimate  of  the 
actual  vorticity  field.  Finally,  this 
method  provides  information  at  only  a 
single  point.  If  information  on  the 
entire  flow  field  is  required,  measure¬ 
ments  must  be  carried  out  sequentially  one 
point  at  a  time.  This  sequential  method, 
although  laborious,  is  straightforward  in 
applications  i,  /o'ving  steady  flows. 
However,  the  method  becomes  very  diffi¬ 
cult,  if  not  impossible  to  implement,  when 
studying  unsteady  flows.  Direct  measure¬ 
ment  of  vorticity  has  been  tried,  for 
instance,  by  injection  of  spherical 
particles  which  rotate  in  the  flow  with  an 
angular  velocity  proportional  to  the  local 
voiticity^  .  Such  methods  suffer  the  same 
drawback  of  insufficient  spatial  resolu¬ 
tion  just  mentioned  and  also  can  be  quite 
complex . 

Recently,  a  novel  velocity  measuiement 
technicjue,  commonly  tefei  •'ed  to  as  Lasei 
Speckle  Velocimetry  (LSV)  or  Paiticle 
Image  Displacement  Velocimetry  (PIDV)  has 
become  available.  This  technique  provides 
the  simultaneous  visualization  of  the  two- 
dimensional  stieamline  pattern  in  unsteady 
flows  as  well  as  the  quantification  of  the 
velocity  field  over  an  entiie  plane.  The 
arivantage  of  this  technique  is  that  the 
velocity  field  can  be  measured  over  an 
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entire  plane  of  the  flow  field  simul¬ 
taneously,  with  accuracy  and  spatial 
resolution.  From  this  the  instantaneous 
vorticity  field  can  be  easily  obtained. 
This  constitutes  a  great  asset  for  the 
study  of  a  variety  of  flows  that  evolve 
stochastically  in  both  space  and  time, 
such  as  the  unsteady  vortical  flows  that 
occur  in  rotorcraft  and  high-single-of- 
attack  aerodynamics.  For  the  background 
of  this  present  technique,  the  reader  is 
referred  to  reference  5. 

The  principle  of  the  technique  is 
given  in  the  ne.xt  section  followed  by  the 
description  of  the  apparatus,  instrumen¬ 
tation  and  procedures.  Section  4  provide 
the  results  and  their  description  for  the 
flow  over  an  high  angle  of  attack  NACA 
0012  airfoil. 


2.  Principle  of  the  Technique 

The  application  of  LSV  or  PIDV  to 
fluid  flow  measurement  involves  several 
steps.  First,  it  is  necessary  to  "create" 
a  selected  plane  or  surface  within  the 
flow  field.  This  is  accomplished  by 
seeding  the  flow  with  small  tracer 
particles,  similarly  to  LDV  applications, 
and  illuminating  it  with  a  sheet  of  coher¬ 
ent  light,  as  shown  in  Figure  1.  A  pulsed 
laser  such  as  a  Ruby  or  a  NdYag  laser,  or 
a  CW  laser  with  a  shutter  is  normally  used 
as  the  light  source.  The  laser  sheet  is 
formed,  for  example,  by  focusing  the  laser 
beam  first  with  a  long  focal  length  spher¬ 
ical  lens,  to  obtain  minimum  thickness, 
and  then  diverging  the  beam  in  one  dimen¬ 
sion  with  a  cylindrical  lens.  The  light 
scattered  by  the  seeding  particles  in  the 
illuminated  plane  provides  a  moving 
pattern.  When  the  seeding  concentration 
is  low,  the  pattern  consists  of  resolved 
diffraction  limited  images  of  the  partic¬ 
les.  When  their  concentration  increases, 
the  images  overlap  and  interfere  to  pro¬ 
duce  a  random  speckle  pattern.  A  multiple 
exposure  photograph  records  this  moving 
pattern.  The  lower  particle  concentration 
originates  a  mode  of  operation  of  the 
technique  referred  to  as  Particle  Image 
Displacement  Velocimetry,  reserving  the 
term  Laser  Speckle  Velocimetry  for  the 
high  particle  concentration  levels  where  a 
random  speckle  pattern  is  actually  formed 
(reference  6).  In  a  second  step  the  local 
fluid  velocity  is  derived  from  the  ratio 
of  the  measured  spacing  between  the  images 
of  the  same  tracer,  or  speckle  grain,  and 
the  time  between  exposures. 

Several  methods  exist  to  con,ett  the 
information  contained  in  the  multiple- 
exposed  photograph,  or  specklegtam,  to 
flow  field  data  such  as  velocity  or  vor¬ 
ticity.  The  recorded  image,  whether 

formed  by  isolated  disks,  in  the  case  of 
low  particle  concentration,  or  speckle 
grains  for  high  particle  concentration  is 
a  complicated  random  pattern.  It  would  be 
very  difficult  to  measure  the  local  dis¬ 
placements  by  visual  or  computer-aided 
inspection.  However,  it  is  important  to 
realize  that  the  multiple  exposure  photo¬ 
graph  results  in  a  periodic  random  image 
from  which  the  periodicity  information  can 


be  retrieved  using  Fourier  or  Auto¬ 
correlation  analysis.  Basically,  the 

multiple-exposed  photographs  or  speckle- 
grams  can  be  analyzed  eitliec  on  a  point- 
by-point  basis,  which  yields  measurements 
of  the  local  displacements  (velocity), 
(refs.  7-fl)  oi  wilh  a  wliole  field  filter¬ 
ing  teclin  i  <pje  ,  wliicli  yields  isovelociry 
contours  (let.  9).  'idle  metliod,  wliicli  lias 
been  selected  and  implemented  by  the  Fluid 
Mechanics  Research  Laboratory  at  the 
Florida  State  University,  is  the  Young's 
fringes  metlioci.  Tlie  local  displacement  is 
determined  using  an  focused  laser  beam  to 
interrogate  a  small  area  of  the  multiple 
exposed  photograph  t ranspa rency .  The  dif¬ 
fraction  produced  by  coherent  illumination 
of  the  multiple  images  in  the  negative 
generates  Young's  fringes,  in  the  Fourier 
plane  of  a  lens,  provided  that  ttie 
particle  images  correlate.  This  is  shown 
schematically  in  Figure  2.  These  fringes 
have  an  orientation  which  is  perpendicular 
to  the  direction  of  the  local  displacement 
and  a  spacing  inversely  proportional  to 
the  displacement.  'The  use  of  Young's 
fringes  eliminates  the  difficulties  of 
finding  the  individual  image  pairs  in  the 
photograph.  The  basis  of  the  Young's 
fringe  method  is  described  in  reference 

a. 

The  photographic  recording  method 
discussed  above  has  the  disadvantage  that 
the  photograph  consists  of  particle  pairs 
v/hich  have  a  180  degrees  ambiguity  in  the 
direction  of  the  velocity  vector.  In 
addition,  it  has  been  shown  (reference  10) 
that  the  velocity  dynamic  range  of  the 
technique  is  limited  to  a  maximum  value  of 
about  10.  In  most  flows  of  interest  (e.g. 
lioundary  layers  and  separated  flows),  this 
dynamic  range  is  not  sufficient  to  capture 
tile  flow  field  in  its  entirety.  These 
limitations  are  critical  when  measuring 
comiilex  flows  liaving  flow  .  reversals  and 
stagnation  areas. 

A  method  to  resolve  both  the  ambiguity 
of  tlie  velocity  vector  as  well  as  to  im¬ 
prove  the  technique's  velocity  range  is 
incorporated  in  this  experiment.  This 
method  proposed  by  Lourenco*'  and 
Adrian'^,  commonly  known  as  "velocity  bias 
technique".  Consists  of  recording  the  flow 
field  in  a  moving  reference  frame,  thus 
superposing  a  known  velocity  bias  to  the 
actual  flow  velocity.  This  effect  may  be 
accompl  i  slied  in  several  ways,  in  particu¬ 
lar,  using  a  moving  camera  during  the 
photographic  recording  or  by  optical  means 
using  scanning  or  rotating  mirrors.  For 
the  data  presented  here,  a  rotating  mirror 
was  used  to  displace  the  image  during  the 
exposure  with  a  p i e-de te rmi ned  velocity. 

A  schematic  of  the  rotating  mirror 
arrangement  is  shown  in  figure  3. 
Consider  two  particle  pairs  A^  and  C 
having  equal  displacements  in  opposite 
directions  in  the  object  plane.  By  intro¬ 
ducing  a  45°  mirror  between  the  camera 
lens  and  the  object  plane,  the  correspond¬ 
ing  displacements  appear  in  the  film  plane 
as  AB  and  CD  with  equal  magnitudes.  When 
the  mirror  is  rotated  by  an  angle  of  ae 
between  exposures,  the  displacements 
corresponding  to  A^ and  D^  appear  in 
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the  film  plane  as  AB‘  and  CD*  with 
different  magnitudes.  The  correct  dis¬ 
placement  or  velocity  with  its  direction 
can  now  be  obtained  upon  removal  of  the 
velocity  bias.  An  example  of  the  flow 
field  obtained  with  and  without  the 
velocity  bias  can  be  seen  in  figure  4. 
This  flow  represent  a  typical  separated 
flow  field  containing  flow  reversal  and 
stagnation  areas  (a  discussion  of  this 
flow  field  is  given  later). 

3.  Apparatus,  Instrumentation  and 
Procedures 


The  time-space  development  of  the 
unsteady  separated  flow  generated  by  an 
high  angle  of  attack  (o=30°)  NACA  0012 
airfoil  impulsively  started  from  rest  is 
examined  using  Particle  Displacement 
Velocimetry.  The  flow  is  created  by 
towing  the  airfoil  in  the  reduced  scale 
Fluid  Mechanics  Research  Laboratory  towing 
tank  facility.  The  tank  is  300  x  200  x 
600  mm.  A  detailed  examination  showed 
that  the  motion  of  the  towing  carriage  is 
smooth  and  vibration  free.  The  towing 
carriage  is  driven  by  a  variable  D.C. 
motor,  and  the  towing  velocity  can  vary 
from  0  to  100  mm/sec.  For  the  photo¬ 
graphy,  a  35mm  camera  (Nikon  F-3)  is  used. 
To  photograph  the  flow  at  regular  time 
intervals,  the  photographic  camera  has  a 
electric  winding  device.  The  photographic 
time  interval  available  with  this  camera 
can  be  continuously  varied  up  to  a  maximum 
of  6  frames  per  second.  Two  options  are 
available  to  fix  the  camera;  one  by 
attaching  it  to  the  towing  carriage, 
which  means  an  observation  point  fixed  in 
relation  to  the  model,  and  the  other  by 
attaching  it  to  the  frame  of  the  water 
tank,  which  means  an  observation  point 
fixed  in  relation  to  the  fluid.  The 
selection  of  these  two  depends  upon  the 
flow  field  being  photographed. 

In  this  experiment  the  airfoil  chord 
is  60ram  and  is  towed  with  a  velocity  of 
22mm/sec.  The  corresponding  Reynolds 
number  was  1400.  The  fluid  used  in  this 
experiment  was  water  seeded  with  4/t/m 
metallic  coated  particles  (ISI  model 
10087).  For  the  illumination,  a  laser 
beam  from  a  5  Watt  Argon-Ion  Laser 
( Spectra-Physics  series  2000)  is  steered 
and  focused  to  a  diameter  of  .3mm  using  an 
inverse  telescope  lens  arrangement.  A 
cylindrical  lens,  with  a  focal  length  of 
-6.34mm,  is  used  to  diverge  the  focused 
beam  in  one  dimensioi',  creating  a  light 
sheet.  The  laser  sheet  is  70mm  wide  and 
illuminates  the  mid-span  section  of  the 
airfoil.  For  the  multiple  exposure,  the 
CW  laser  beam  is  modulated  using  a  Bragg 
cell.  In  this  experiment,  the  laser  power 
density,  I^,  ,  of  the  sheet  was  .27  W/mm^  . 
In  order  to  record  the  time  development  of 
the  flow  field,  the  camera  was  attached  to 
the  towing  carriage  and  the  frequency  of 
which  the  multiple  exposures  were  taken 
was  set  at  2.0Hz.  The  aperture  of  the 
lens  with  a  focal  length  of  50mm  and  a 
spacer  of  12mm,  was  set  at  Ftt5.6  and  the 
resulting  magnification  factor  was  0.40. 
For  the  multiple  exposure,  the  time 
between  exposures,  T,  and  the  exposure 


'time,  t,  are  chosen  according  to  the  cri¬ 
teria  discussed  in  reference  6.  The  time 
between  exposures  was  10msec.  For  optimum 
exposure,  the  exposure  time  was  1msec, 
whicli  corresponded  roughly  to  (d,/MVmax), 
where  D  is  the  analyzing  beam  diameter,  V 
is  the  maximum  expected  velocity  plus  the 
shift  velocity  in  the  field  and  d  is  the 
paiticle  image  diameter  expressed  in  terms 
of 

1^ 

d  =  ( d^  +  d^  )  ^ 

with  d  the  particle  diameter  and  d^  the 
edge  s^iread  caused  by  the  limited  response 
of  the  recording  optics  (ref.  6). 

Data  Processing 

The  fringe  images  were  acquired  and 
analyzed  using  tlie  digital  image  analysis 
system  of  tlie  Florida  State  University 
FMRL  (Fig.  5).  This  system  consists  of 
the  following  components:  a  DEC  LSI-11/73 
host  computer,  Gould  IP-8500  Digital  image 
processor  which  includes  four  memory  tiles 
for  storage  of  image  data  in  a  512  x  512 
format  with  a  resolution  of  8  bit  per 
pixel,  a  frame  digitizer,  a  pipeline  pro¬ 
cessor  and  a  video  output  controller  to 
convert  digital  to  analog  information  for 
display  on  a  color  monitor.  The  system 
also  includes  a  two-dimensional  Klinger 
traversing  mechanism  with  a  controller  for 
tlie  purpose  of  automatically  scanning  the 
film  transparencies.  Two  methods  are 
available  and  used  for  fringe  analysis 
(ref.  10).  The  first  one  is  an  inter¬ 
active  method  in  the  sense  that  it 
requires  the  assistance  of  an  operator. 

The  inconvenience  of  the  one¬ 
dimensional  averaging  method  is  the  need 
for  an  external  adjustment  of  the  angle  of 
the  fringes  by  an  operator.  This  problem 
can  be  by  passed  by  computing  the  velocity 
components  along  independent  directions. 
Because  each  line  of  the  fringe  frame  can 
be  considered  as  a  noisy  periodic  signal 
with  variable  phase,  the  automatic  deter¬ 
mination  of  a  velocity  component  can  be 
performed  only  by  averaging  over  a  quanti¬ 
ty  independent  of  the  phase.  The  autocor¬ 
relation  for  each  line  or  its  Fourier 
transform  for  the  power  spectrum  satisfies 
t)iis  requirement.  The  m  velocity  compo¬ 
nent  can  be  computed  from: 


g(  u ) 


511 

L 

n=0 


[  I  (  m,  n  )  I  (  m+u  ,  n  )  ) 
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-511<u<511 


This  algorithm  has  been  implemented 
using  the  pipeline  processor  of  the  Gould 
IP-8500  image  processor  to  perform  simul¬ 
taneously  the  autocorrelation  for  all  the 
lines  of  a  frame.  For  an  accurate  esti¬ 
mate  of  the  velocity  magnitude  and  direc¬ 
tions,  four  of  such  full  image  operations, 
yielding  four  autocorrelation  functions, 
are  required.  From  these  the  velocity 
vector  is  determined  by  selecting  the 
values  of  the  components  which  have  been 
comouted  from  autocorrelations  having  the 
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highest  SNR,  and  visibility.  (The  compu¬ 
tation,  which  includes  the  determination 
of  the  fringe  angle  and  position  updating 
of  the  film  transparency  scanning  mecha¬ 
nism,  is  completed  in  a  few  seconds, 
typically  4-5  sec,  using  the  PDF  11-73 
compute  t . ) 

The  overall  accuracy  of  the  technique 
was  evaluated  using  a  method  described  in 
reference  10.  A  uniform  flow  field  is 
created  by  producing  a  multiple  exposure 
photograph  of  the  still  seeded  water,  in 
the  water  tank,  with  a  camera  moving  at 
constant  speed.  For  the  multiple  exposure 
photograph  a  number  of  time  between 
exposures  are  used,  thus  resulting  in 
photographs  with  particle  pairs  at 
different  known  distances. 

In  the  absence  of  a  systematic  bias, 
the  standard  deviation  of  the  obtained 
velocity  distribution  is  an  estimate  of 
the  mean  measurement  error.  From  the 
error  analysis,  it  is  believed  that  the 
velocity  data  is  obtained  with  an  accuracy 
of  2  percent  or  better. 

4.  Results  and  Discussion 

Typical  multiple  exposure  photographs 
of  the  flow  generated  by  the  impulsively 
started  NACA  0012  from  rest  for  different 
times  are  shown  in  figure  6.  The  photo¬ 
graphic  arrangement  was  purposely  adjusted 
to  enhance  the  view  of  the  flow  field  on 
the  upper  surface  of  the  airfoil  rather 
than  to  show  the  entire  flow  around  the 
airfoil.  Consequently,  the  details  of  the 
flow  under  the  airfoil  can  not  be  seen 
clearly  in  these  photographs.  The  angle 
of  attack  of  the  airfoil  is  set  at  30°. 
These  pictures  display  the  flow  field  from 
the  leading  edge  to  a  downstream  location 
of  about  1.5  chords.  Photographs  were 
taken  at  a  frequency  of  2Hz.  A  total  of 
34  pictures  were  obtained  covering  the 
range  of  t‘  from  0.1  to  5.  The  non- 
dimensional  time,  t'  (»  Ut/c,  where  U  is 
the  free  stream  velocity,  t  is  the  time 
between  two  successive  pictures,  and  c  is 
the  chord  of  the  airfoil)  between  suc¬ 
cessive  pictures  was  0.167.  However,  in 
figure  6,  only  a  limited  number  of  pic¬ 
tures  are  included.  The  quadruple  exposed 
photographs  shown  here  increase  the  SNR 
(signal  to  noise  ratio),  the  fringe  visi¬ 
bility,  and  provide  an  excellent  flow 
visualization.  The  added  advantage  of 
providing  a  good  flow  visualization  is  an 
asset  of  the  PIDV  technique. 

From  the  flow  visualization  pictures 
shown  in  figure  6,  the  following  obser¬ 
vations  are  made.  At  the  start  of  the 
airfoil  a  vortex,  at  the  trailing  edge, 
commonly  known  as  "starting  vortex",  is 
generated  and  is  carried  away  from  the 
body.  concomitant  with  this  is  the 
generation  of  a  separation  bubble  at  the 
leading  edge  of  the  airfoil  (figure  4a). 
At  a  later  time,  for  example  at  t‘  -  1.2 
(figure  4b),  the  separation  bubble  grows 
into  an  isolated  primary  vortex  with 
"secondary  vortices"  following  behind  it. 
Similar  type  of  vortex  structure  was  also 
observed‘“in  the  flow  behind  a  circular 


cylinder.  This  multiple  vortex  structure 
continue  to  grow  together  until  the  t‘ 
reaches  a  value  of  about  2.5  (figure  4c 
and  4d).  At  t'  =  2.5  (figure  4d),  because 
of  the  close  proximity  of  the  primacy 
vortex,  a  trailing  edge  vortex  is 
generated.  At  t‘  =  2.75  die  primary 
vortex  afjruptly  moves  away  from  tin;  sur¬ 
face  of  the  airfoil  leaving  behind  a 
"vortex  sheet"  like  structure  (figure  4e ) . 
For  t‘  >  3.0,  this  "vortex  sheet"  rolls  up 
into  distinct  vortices  and  they  grow  in 
size  with  time  as  shown  in  figure  4f  -  4q. 
During  this  process,  the  trailing  edge 
vortex  also  grows  and  as  a  result  the 
whole  flow  field  becomes  very  complex. 
Close  to  the  surface  of  the  airfoil,  a 
small  vortex  remain  present  for  t’  >  3.0. 
This  vortex  has  the  same  sign  of  rotation 
as  ■.  -■  trailing  edge  vortex.  A  similar 
vorteA  structure  was  observed  by  Ho' \  who 
calls  it  an  "induced  vortex"  and  associ¬ 
ates  it  with  unsteady  separation 
phenomenon . 

The  velocity  data  is  acquired  in  a 
Cartesian  mesh  by  digital  processing  of 
the  Young's  fringes,  produced  by  point-by¬ 
point  scanning  of  the  positive  contact 
copy  of  the  photograph.  The  scanning  step 
size  and  the  dimension  of  the  analyzing 
beam  are  0.5mm,  which  corresponds  to  a 
spatial  resolution  of  about  1.25mm  in  the 
object  plane  or  about  0.02c.  The  fringes 
were  processed  using  the  method  described 
in  the  previous  section.  The  resultant 
two-dimensional  velocity  fields,  corre¬ 
sponding  to  figure  6a  -  6h,  are  shown  in 
figure  7.  The  length  of  each  vector  is 
proportional  to  the  local  velocity  at  that 
point.  The  color  code  superimposed  on  the 
velocity  data  represents  the  vorticity 
level,  the  magnitude  of  which  is  given  by 
the  color  bar  at  the  bottom  of  the 
picture.  The  red  and  green  colors  repre¬ 
sent  the  peak  positive  and  negative  vor¬ 
ticity  regions  respectively.  This  type  of 
display  clearly  depict  the  various  regions 
of  vorticity  and  its  strength.  As  dis¬ 
cussed  above,  the  presence  of  primary 
vortex,  secondary  vortices,  vortex  sheet 
trailing  edge  vortex  and  the  induced 
vortex  on  the  surface  of  the  airfoil  are 
clearly  depicted  along  with  their  time- 
space  development.  The  detailed  analysis 
of  this  data  is  being  conducted  at  this 
time  and  will  be  reported  later. 

5 .  Cone lus i ons 

A  recently  developed  velocity  measure¬ 
ment  technique,  known  as  Particle  Image 
Displacement  Velocimetry  (PIDV),  )ias  been 
briefly  described.  Using  this  technique, 
the  time-space  evolution  of  the  flow 
generated  by  an  impulsively  started  high 
angle  of  attack  (a=30°)  NACA  0012  airfoil 
is  presented.  This  experiment  illustrated 
the  technique's  capabilities  to  record 
with  accuracy  the  complex  unsteady 
separated  flow. 

The  technique  has  been  shown  to  pro¬ 
vide  both  flow  visualization  and  quanti¬ 
tative  measurements,  which  include  the 
velocity  and  vorticity  fields. 
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The  unsteady  separated  flow  field 

generated  by  an  high  angle  of  attack  air¬ 
foil  contains  many  large  scale  vortical 
structures  such  as;  primary  vortex 

generated  at  the  leading  edge  with  second¬ 
ary  vortices  upstream  of  it,  trailing  edge 
vortex,  vortex  sheet  and  an  induced  vortex 
in  the  upper  surface  of  the  airfoil.  The 
origins  and  time  development  of  tliese  ate 
clearly  depicted  by  the  instantaneous 

velocity  and  vorticity  fields  obtained 
using  PIDV. 
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Figure  1,  Schematic  arrangement  for  the  photographic  recording. 


Figure  2.  Schematic  arrangement  for  obtaining  Young's  fringes. 
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Figure  3.  Rotating  Mirror  arrangment  for  the  Velocity  Bias  Technique 
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Figure  4.  Instantaneous  Velocity  Field  a)  Before  removal  of  velocity 
bias  b)  after  removal  of  the  velocity  bias. 


Figure  5.  Schematic  of  the  data  analysis  system 


neous  velocity  and  vorticity.  fields;  a)  J  -0.3;  b) 


APPENDIX  n 

Fast,  Adaptive  Summation  of  Point  Forces  in  the  Two-Dimensional  Poisson  Equation, 
Journal  of  Computational  Physics  Paper. 


I. 


JoLirnal  of  CompuUiiional  Physics  -  4017 


JOLKNAL  ()l  ( OMIMJTATIONM.  I'H'i'ilCS  83.  18X)-0(X)  (1989) 


Fast,  Adaptive  Summation  of  Point  Forces  in 
the  Two-Dimensional  Poisson  Equation 

Li:ON  VAX  DdMMI'.LP.N 
hepat  tnu'n!  of  Mcduitticol  F.n,i;ith‘crin\;. 

f  A  \f  (i  /  FSU  CnlltXt’  f>l  l.n^inrcrinii  <<<  Sitficrc<>tHfUilcr  Cff/npuldHom  lir^canh  hi'Hiiuir. 
/  ioni/,t  Slatr  ( 'nircrsiiy^  TtiUaluisscc.  Flortiia 


AND 


[  j  Kl  A.  RUNDI-.N.STrilNCR 

Ocptirrmrnl  nf  Cotupuft-r  Science, 
ilnrut.i  Slate  i'tttvcr.stiy.  Tailnhaisee.  Florida 

lU'cci'.ccJ  June  2,  1987;  rcvi.Accl  .September  .''<1  I9SS 


Direct  .summ.ilioii  of  tlio  'clocily  Held  jiilrotiiiecd  by  point  vorlice.'  lends  lo  be  lime 
eonsummg  since  the  'cloei:\  of  each  vorte.x  is  found  as  a  sum  over  all  other  vortices.  The 
rc.sultiiiu  miinbcr  of  mimen^.il  operations  is  proponional  to  the  square  of  the  number  of 
lorticcs  Here  a  relatively  Miiipic  prticediirc  is  outlined  uliieh  significanlly  reduces  the  number 
of  operalions  by  rc|il.ieinc  sciecied  partial  sums  by  asymptotic  scries.  The  resulting  number  of 
operalions  appears  tr'  vary  roughly  in  proportion  to  the  number  of  unknowns,  corresponding 
to  a  "f,l.^i"  solver.  •  I'JSV  I'tcw.  Inc 


I.  Intkoouciion 

liicoinpfcs.sible  flow  tit  high  Reynolds  number  with  Kirgc-sctiic  separtuion  can  be 
difllciill  III  compute  .since  ihe  vorlicity  lends  to  concentrate  m  limiied  parts  of  the 
flow  Held.  Vortex  methods  [I]  attempt  to  reduce  the  number  of  variables  needed 
for  the  computation  by  describing  only  the  vorticity,  in  its  simplest  form,  by  a  senes 
of  delta-functions  or  point  vortices: 


n; 


=  X  /Vdlx-x,).  (1) 


The  flow  velocity  is  related  to  the  vorticity  by  the  solu^.on  of  a  Poisson  equation, 
wilh  the  vorticity  as  forcing  function,  resulting  in  the  stream  function.  The  flow 
velocity  is  found  by  t.ikmg  the  curl  of  the  strctim  function.  The  solution  for  ;i  senes 
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of  delta  functions  can  be  found  aiul  leads  to  the  following  expression  for  the  flow 
velocity: 

=  i='.2 . V,  (2) 

1=1  ^1  ^ : 

I  •  / 


where  Z  is  the  complex  position  .v  +  —  1  r  and  K*  the  complex  conjugate  \elocity 

1/  -  ^  —  I  r.  In  general,  it  will  be  necessary  is'  add  to  this  flow  velocity  a  solution 
of  a  Laplace  problem  to  take  care  of  the  boundary  conditions. 

While  ihe  sum  in  Eq.  l2)  may  easily  be  evaluated,  the  number  of  terms  is  propor¬ 
tional  to  the  square  of  the  number  of  vortices  /V.  Thus,  the  computational  effort 
increases  rapidly  when  the  number  of  vortices  increases.  In  contrast,  various  mesh- 
b.iscd  solution  procedures  for  the  Poisson  equation  arc  able  to  find  the  solution  in 
a  computational  time  roughly  proportional  to  the  number  of  mesh  cells.  As  a  result, 
a  point  \ortex  descri|Hion  seems  most  useful  if  (a)  the  number  of  vortices  is  much 
smaller  th;m  the  number  of  mesh  cells  needed  to  describe  the  flow  (i.e.,  the  vorticity 
is  restricted  to  a  small  part  of  the  total  domain);  (b)  the  point-singularity 
description  itself  is  of  p.irticular  interest  and  the  errors  induced  by  a  mesh-based 
rei'resenlation  must  be  avoided;  or  (c)  the  infinite  domain  implicit  in  Eq.  (2)  is  to 
be  preserved.  Certainly  discrete  vortex  representations  have  drawn  and  continue 
to  draw  considerable  theoretical  and  numerical  interest.  In  addition,  the  Poisson 
equation  is  not  unitiue  to  fluid  mechanics;  it  arises  in  other  fields  such  as  electro¬ 
magnetism  and  gravitation.  For  these  reasons,  more  efficient  procedures  to  evaluate 
the  solution  under  point's isc  forcing  arc  of  considerable  interest. 

Various  methods  to  reduce  the  computational  effort  h.ive  been  proposed. 
Anderson  [2]  used  ;i  last  Fourier  transform  mcthotl.  with  corrections  for  the 
interactions  between  nearby  vortices.  However,  some  of  the  mentioned  advantages 
of  the  vortex  mcthoil  are  lost  due  to  the  presence  of  the  mesh  In  addition,  for  high 
accuracy  the  evaluation  of  the  interactions  between  neighboring  vortices  can 
become  computationally  intensive. 

,‘\n  alternative  approach  followed  in  this  paper  is,  to  group  the  vortices  sptiiially 
and  to  api'roximate  the  elTccts  induced  by  each  group  at  larg  listances.  Appel  [3] 
and  Barnes  and  Hut  [4]  made  approximations  using  a  single  rcphiccmenl  element. 
Yet,  using  such  approximations,  high  accuracy  is  difficult  to  achieve  while  the 
algorithm  tends  to  be  scalar. 

In  contrast,  the  present  study  uses  a  L.iurent  series  ai'proximation  for  the 
velocity  induced  by  ctich  group.  This  approximation  t.ikcs  the  form 
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where  Z„  IS  u  svi\t.ibly  chosen  origin  point  ior  the  group  oi  vortices  ;ind  the  sum 
m  (3b)  extends  over  all  vortices  in  the  considered  group.  The  Laurent  series  allow 
the  desired  accurticy  to  be  maintained  by  the  choiee  of  the  truncation  of  the  infinite 
sum.  [n  iiddition,  when  the  point  /  at  which  the  velocity  is  to  be  evaluated  is 
sufTicienily  far  distant  from  the  group  of  vortices,  the  series  converges  geometrically 
and  only  a  limited  number  of  terms  is  needed  for  given  accuracy.  Savings /?  com¬ 
putational  effort  result  when  the  number  of  terms  nccticd  for  the  Laurent  series  is 
sufficiently  small  eomptued  to  the  number  of  vortices  m  the  group.  For  that  reason, 
a  minimum  group  size  exists  below  which  further  savings  are  not  made.  Using  the 
adaptive  algorithm  on  a  CYBER  205  computer.  Van  Dommelen  and  Rundensteiner 
[5]  found  that  this  group  size  is  of  the  order  of  a  100  elements. 

At  about  the  same  time,  similar  ideas  were  developed  by  Rokhlin  [6]  and 
(ircengaid  and  Rokhlin  [7]  In  fact,  an  adaptive  algorithm  developed  by  Carrier. 
Circengard,  and  Rokhlin  [X]  is  quite  similar  to  the  present  one  in  both  the  use  of 
Laurent  series  and  the  grouping  involved.  An  important  difference  between  the 
procedures  is  how  the  adaptive  group  structure  is  addressed.  While  the  procedure 
[S]  is  based  on  five  topological  sets  expressing  the  relationships  between  groups, 
the  present  procedure  is  based  on  an  unusual  numbering  system  of  the  groups.  The 
numbering  system  is  generated  simuitancously  with  the  group  structure;  it  leads  to 
a  relatively  simple  and  streamlined  program  logic. 

The  procedure  of  Grccngard  and  Rokhlin  recasts  the  Laurent  scries  as  Taylor 
series  to  achieve  further  reductions  in  computational  operations,  an  enhancement 
not  yet  incorporated  in  the  present  scheme,  liowc'  t,  unless  the  number  of  vortices 
is  sufficiently  large,  the  possible  savings  seem  limited.  Furthermore,  not  recasting 
the  scries  offers  some  compensating  advantages,  such  as  reduced  storage  (only  a 
s .vnishinjily  sivudi  fraction  of  the  Laurent  senes  exp.insions  need  be  stored), 
mcreascjf’  vector  length,  and  less  overhetid. 

In  Its  present  form,  our  proeevlurc  can  he  divided  into  two  p;irts:  generation  of 
an  adaptive  panel  structure,  to  groups  the  vortices  sp.iiially,  ami  determination  of 
the  velocity.  The  next  two  sections  viescribc  each  of  these  steps  m  turn. 

2  (jIM  KAIION  AND  Ni'MIII  ItlN't;  OI  titl  I’.ANI  IS 
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In  order  to  use  Laurent  scries  effectively,  the  vortices  must  be  sp.ititilly  groupevt 
together  Figure  I  illustrates  a  typical  grouping  for  the  ex.imple  of  flow  about  a 
circular  cylinder.  In  this  example,  there  are  16.479  vortices  outside  the  evliiulci 
(shown  as  dots)  ;iikI  an  equal  number  of  mirror  vortices  iii'-nle  the  cyliiKlei  I  not 
shown ). 

The  procedure  for  generating  this  panel  sliuettirc  is  shown  in  big.  2  fhe  first  few 
steps  are  further  illustrated  m  Fig.  3.  The  starling  domain  is  taken  as  the  smallest 
square  that  encloses  all  vortices.  This  sv|uaie  is  subvlivideri  into  four  squares,  or 
subpancls,  of  equal  size  (indicated  as  A.  B.  CL  and  1)  in  big  3).  The  vortices  arc 
reordered  so  as  to  group  the  vortices  in  each  of  the  four  subpanels  together  (iti  the  ^ 


f'lG.  I.  lixainplc  I'.iiicl  slrav'Hifc  gcacralcd  for  flow  aboiil  a  toMilar  cvlimlcr  The  I6,(XK1  xorliccs 
cHilside  the  cylinder  arc  shown  as  dot.,  an  C(|ual  nundicr  of  mirror  -.orlices  wiihin  the  cylinder  are  not 
show  n. 


C’YBLR  205  impIcmciUiilion,  the  bullt-iii  vcclor  funciuMi  QSX’CMPRS  was  used  for 
this  pur|5ose). 

Information  iibout  ctich  of  the  panels  is  stored  m  memory.  This  mforimition 
includes  (he  position  of  the  panel  and  the  storage  loeaiions  of  the  llrst  and  the  hist 
vortex  within  the  panel  It  ;dso  includes  an  identifying  ptincl  number  dcTincd  later. 

After  subdivision,  execution  transfers  to  the  first  of  the  four  subpanels  genertited, 
;ind  a  decision  is  imide  vvhether  this  ptincl  should  be  further  subdivided  The 
decision  is  based  on  the  number  of  vortices  in  the  panel;  if  sufficient  vortices  arc 
present,  for  example,  more  than  lOO,  the  panel  is  further  subdivided,  Metinwhile  the 
panel  datct.for  the  remtiming  three  subpancls  is  temporarily  stored  awtiy  m  a  last-m, 
fir|^-out  buffer,  (The  last  ptiit  of  the  memory  allocated  for  the  panel  information 
was  used  as  buffer.) 
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I  ic.,  2  D.isic  flow  chan  for  generating  the  p.uKl  slmclurc. 


If  the  eonsidcreri  panel  docs  not  eoiitain  stifncicnl  '  ('rliees.  the  next  panel  will  be 
retrieved  from  the  buffer  fo;  possible  fiiithcr  subdivisions.  Proceeding  in  this 
manner,  the  entire  domain  is  subdivided  into  panels  eontaining  a  limited  number 
of  vortiees. 

Each  of  the  panels  ia  given  a  unique  number  to  simphh’  identification  of  the 
panel  and  its  place  in  the  panel  structure.  Figure  3  illustrates  thtit  the  storage  is 
always  kept  in  order  of  increasing  panel  number.  The  actual  definition  of  the  panel 
number  is  illustrated  in  I'lg.  4:  all  panels  which  could  be  created  by  the  subdivision 
process  can  be  represented  as  a  series  of  uniform  divisions  of  the  domain  For  each 
I'f  these  uniform  subdivisions,  the  .v-  and  r-positions  can  be  given  a  binary  number. 
The  four  panels  generated  at  the  first  level  of  subdivision  ean  be  numbered  using 
a  onc-digit  binary  number  (top  of  Fig.  4).  Each  additional  level  of  subdivision 
requires  one  additional  digit. 

Therefore,  the  binary  digits  determine  the  position  of  the  panel  The  number  of 
binary  digits  determines  the  subdivision  level.  It  follows  that  the  binary  digits  of  the 
.V-  and  1 -positions  dcsciilte  the  panels  uniquely.  The  complete  information  is  stored 
111  a  single  panel  numbei  using  the  following  proeediirc:  ij' 


(a)  Increment  each  digit  in  holli  the  .v-  and  v-position  by  one.  so  that  binaiv 
zero  becomes  1  and  binary  one  becomes  2. 

(b)  “Interletive”  the  resulting  digits  of  the  .v-  and  t-pnsitions  into  a  single 
number,  so  that  the  odd  digits  become  the  digits  of  the  .v-position  and  the  even 
digits  those  of  the  r-position. 

(e)  Add  trailing  zeros  to  obtain  a  final  panel  number  with  a  fixed  and 
predetermined  number  of  digits.  The  205  procedure  cliooscs  a  28  tiigit  panel 
number. 

The  procedure  is  illustrated  in  Fig.  4  for  example  panels.  Since  the  highest  value 
of  the  digits  in  tire  obtained  panel  number  is  2,  it  can  be  considered  as  the  represen¬ 
tation  of  a  number  m  a  ha.se-3  notation. 
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I-IG,  3.  Mic  first  tlircc  stepv  m  gcncraltng  the  panel  structure  of  i'lg  !  Slioun  arc:  ihc  subsequent 
(Jnisions  of  liic  tIomaiM.  the  of*lcr  m  which  the  vortices  arc  stored.  :lic  onlet  m  which  the  infornialion 
about  the  panels  is  stored,  and  the  nuinbcnng  of  the  panels  ^ 
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I  If;.  4  Dertnilion  of  llic  panel  ntiinbci  of  example  |iaiiclx 

By  construction,  tlic  panel  minibcr  contains  all  the  infornialion  ahoul  the  panel: 
the  non-zero  digits  determine  the  panel  position  and  the  number  of  pairs  of  non¬ 
zero  digits  the  subdivision  level.  I’articularly  important  pro|K'Hies  arc: 

(i)  For  any  given  panel,  the  panel  numbers  of  neighboring  panels  (if  the 
same  size  may  be  round  by  simply  binary  manipulations,  (l-oi  e.xample,  to  ntul  (he 
panel  at  the  same  r-|iosuion  but  the  previous  .v-position,  dcicrmme  the  odd  non¬ 
zero  digits  of  the  panel  number,  giving  the  v-location.  and  do  a  binary  sutilracdon 
of  unity  to  find  the  digits  of  the  sough;  |iancl  number  1 


I  3  0  — ►  1 

y  =  1  -►  1 

II 

panel  number:  1  1  0  0  0  0  0  0...0  0 


(ii)  For  any  given  panel  number,  tlie  panel  number  of  tlie  next  larger 
“motlier”  panel  containing  the  given  panel  is  found  by  setting  tlie  last  two  non-zero 
digits  to  zero,  (The  last  two  non-zero  digits  were  generated  by  the  last  panel 
subdivision.) 

(iii)  Subpanels  of  any  panel  have  a  panel  number  greater  tha:’  the  original 
panel,  but  less  than  the  next  panel  of  equal  size.  (The  subpanels  have  the  same 
leading  digits  as  the  original  panel,  but  non-zero  trailing  digits.)  This  property 
implies  that  in  order  of  increasing  panel  number,  panels  are  arranged  in  “families," 
with  the  subpanels  always  immediately  following  the  panels  of  which  they  arc  a 
part. 

The  algorithm  for  generating  the  panels  described  at  the  start  of  this  section 
generates  them  in  order  of  increasing  panel  number,  subdividing  the  current  lowest 
panel  before  moving  on  to  the  next  panel. 

Since  each  subdivision  adds  two  more  non-zero  digits,  the  total  number  of  digits 
in  the  panel  number  limits  the  smallest  panel  that  can  be  defined.  In  the  205 
implementation,  this  total  number  of  digits  was  chosen  to  be  28,  since  28-digit 
numbers  arc  the  largest  basc-3  numbers  than  can  be  stored  in  ;i  single  205  memory 
location,  saving  storage  and  computational  operations.  In  28  digit  representation, 
the  smallest  ptincl  c<m  be  about  16,000  times  smaller  thtin  the  original  domain, 
which  would  seem  sufficient  for  most  purposes. 


.V  DftTiatMINATION  Ol'  THI-.  Via.DCITY 

The  velocity  is  (.letcrmincd  in  a  single  pass  over  till  panels  in  the  order  in  which 
they  were  gcncratc(.l  as  described  in  the  previous  section,  fhe  procedure  is  outlined 
m  Fig.  5 

For  each  panel,  a  "neighbiirhood"  of  vortices  is  cstablishcil,  consisting  of  the 
vortices  both  within  the  panel  itself  and  in  the  ptincls,  of  at  least  equal  size,  sharing 
a  boundary  line  or  a  corner  point  with  the  considered  panel  The  vortices  in  this 
neighborhood  arc  not  summed  by  the  Laurent  scries  expansion  of  the  considered 
panel.  This  restriction  ensures  that  the  Laurent  scries  converges  exponentially. 
Instead,  in  evaluating  the  velocity  induced  on  the  neighborhood,  the  original  sum 
in  Ec|.  (2)  is  used.  This  sum  is  only  performed  for  panels  which  arc  not  further 
subdivided;  for  panels  which  arc  subdivided,  the  velocity  is  evaluated  by  means  of 
the  subpanels. 

Laurent  scries  can  be  used  for  all  vortices  outside  the  neighborhood  of  a  panel. 
However,  to  reduce  the  computational  effort,  the  Laurent  series  is  only  used  for 
those  vortices  which  cannot  be  evaluated  by  means  of  the  Laurent  series  of  the  next 
larger  "mother”  panel:  the  single  Laurent  series  of  the  mother  panel  is  more 
efficient  than  the  four  Laurent  scries  of  its  subpanels.  Therefore,  the  Laurent  senes 
of  any  panel  is  used  only  for  the  vortices  within  the  neighborhood  of  the  mother, 
but  outside  the  neighborhood  of  the  panel  itself 


5.  Il.isic  flow  cliarl  for  ihc  iJclcrminalion  of  Ihc  vclocily. 

In  tins  scheme,  al  each  stage  tlie  smallest  possible  number  of  Laurent  series  is 
used,  for  /V  vortices  resulting  in  the  C>(/Vin  A')  operation  counts  of  the  ne.xt  section. 
In  the  jirocedure  of  Grcengard  and  Rokhlin  [7]  this  operation  count  is  further 
reduced  to  0(N)  by  recasting  the  Laurent  series  as  Taylor  .series;  however,  the 
present  procedure  has  the  advantage  of  being  less  complex  and  requires  only  a 
single  sweep  over  the  panel  structure  to  evaluate  the  velocity. 

To  incorporate  tlie  Taylor  series  within  the  present  procedure,  the  evaluation  of 
the  neighborhood  of  the  mother  would  have  to  mollified.  For  each  suitable  panel 
within  this  neighborhood,  the  sum  (3a)  would  be  replaced  by  a  recasting  of  the 
Laurent  series  into  a  Tailor  .series.  Additional  steps  would  be  needed  to  transfer  the 
Taylor  scries  of  the  larger  panels  to  the  subiianels  and  to  tidd  the  contributions  of 
these  series  to  the  velocity. 

Clearly,  this  will  increase  program  complexity  and  scalar  overhead.  In  addition. 
It  requires  that  the  neighborhood  of  the  mother  is  described  in  terms  of  individual 
panels.  I'he  present  procedure  describes  this  neighborhood  in  terms  of  a  small 
number  of  vectors  of  vortices,  increasing  vcctorization.  Furthermore,  the  present 


procedure  has  a  slorage  advantage:  wlicn  the  final  subpancl  of  any  group  of  4  is 
reached,  the  four  I.aurciU  series  of  lire  subpaircls  can  he  combined  into  the  l.auieiU 
series  of  the  motlicr  (bottom  of  Fig.  5).  Tlic  four  Laurent  senes  of  the  subpanels  can 
then  be  discarded;  they  arc  no  longer  needed.  As  a  result,  at  any  time  only  a  small 
fraction  of  the  Laurent  scries  need  be  stored.  On  the  other  hand,  using  Taylor 
series,  no  obvious  way  to  avoid  storing  the  Taylor  series  coefficients  for  each  panel 
is  evident.  This  c:in  be  a  disadvantage  since  each  series  represents  a  set  of  coef¬ 
ficients  while,  in  addition,  the  total  number  of  panels  may  be  difficult  to  estimate 
precisely  beforehand. 

In  the  actual  implementation  of  the  procedure  in  Fig.  5,  the  first  step  is  identifica¬ 
tion  of  the  neighborhood  of  each  panel.  The  present  procedure  sttirts  out  by 
identifying  the  individual  digits  of  the  binary  .v-  and  r-positions  of  the  panel.  By 
performing  unit  binary  additions  and  subtractions,  the  panel  numbers  of  the  eight 
neighboring  panels  of  the  same  size  arc  found.  For  each  of  these  eight  panel 
numbers,  the  corresponding  panel  is  located.  In  case  any  of  the  eight  panels  is 
undefined,  the  panel  with  the  largest  panel  number  less  than  or  equal  to  the  sought 
one  is  selected.  On  behalf  of  the  properties  of  the  panel  number,  the  selected  panel 
will  always  enclose  the  sought  panel,  ensuring  the  geometric  convergence  of  the 
Laurent  series.  Since  the  panel  numbers  arc  ordered,  an  appropriate  search  on  a 
scalar  machine  is  binary;  the  CYBER  205  implementation  switches  to  the  vector 
function  Q8SLT  when  the  search  interval  c.xicnds  over  less  than  500  panels. 

After  the  neighboring  panels  have  been  located,  the  storage  locations  of  the  vor¬ 
tices  m  the  neighborhood  arc  simply  the  combination  of  the  storage  locations  of  the 
vortices  in  each  of  the  nine  panels.  The  subdivision  process  of  the  previous  section 
reordered  the  vortices  so  as  to  group  vortices  in  the  same  panel  in  contiguous 
storage  locations  or  vectors.  As  a  result,  the  neighborhood  is  described  by  at  most 
nine  vectors  of  vortices,  and  an  additional  check  is  made  to  identify  contiguous 
sectors  which  can  be  described  by  a  .single  vector.  (In  particular,  the  four  subpanels 
of  the  larger  panel  containing  the  considered  panel  describe  a  single  vector  of 
lortices  )  Since  the  number  of  vortices  per  panel  is  never  small,  the  computations 
remain  efficient  on  the  205. 

The  next  step  in  the  procedure  in  Fig.  5  is  the  evaluation  of  the  original  sum  in 
llq.  (2)  for  panels  which  are  not  further  subdivided.  This  sum  was  split  into  real  and 
imaginary  parts  and  modified  to: 
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These  c.vprcssions  ai 

re  equivalent 

to  the  original  sum  m  Lc],  (2)  when  the  ' 

value  of 

(i,  vanislics.  They  arc  cc|uivalcnt  to  tiic  sum  m  Eq.  (2)  to  macliinc  precision  when 
the  value  of  r/,  equals  tlic  machine  epsilon.  The  addition  of  the  f/,-tcrm  in  (4a)  and 
(4b)  has  the  advantage  of  avoiding  the  singularity  in  the  /  =  /  term  while  limiting 
the  effect  of  numerical  inaccuracy.  For  larger  values  of  il,,  the  velocity  corresponds 
to  vortices  with  nnitc  core,  which  tend  to  improve  the  numerical  properties  of  a 
vortex  representation  [9].  Expressions  more  elaborate  than  Eqs.  (4a)  through  (4c) 
could  be  used  [10];  since  they  increase  the  computational  time  for  the  original 
algorithm,  they  are  likely  to  enhance  the  relative  performance  of  the  present 
algorithm. 

The  coefTicients  of  the  Laurent  scries  follow  from  the  sum  (3b).  The  eocfTieicnts 
may  be  split  into  real  ;ind  imaginary  parts  .•)*.  and  Ictuiing  to  the  following 
recursive  relationships: 
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(5e),  (5f) 

In  order  to  avoid  possible  inaccuracy  caused  by  underflow  of  terms,  the  .v-  and 
y-positions  were  measured  from  the  center  of  the  panel  and  normalized  with  half 
the  linear  panel  dimension. 

For  the  evaluation  of  the  neighborhood  of  the  motltcr  panels  in  Fig.  5.  the 
Laurent  series  {3a)  is  used,  ^it  into  real  and  imaginary  parts,  the  scries  can  be 
written:  ' 
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It  can  readily  be  sliown  that  the  terms  induced  by  any  individual  vortex  i 
converge  geometrically  with  a  convergence  ratio 


|t»l  _  2  Z() 

It*  (.  1 1  ~  Zf) 


(7) 


if  Z,|  is  the  position  of  the  center  of  the  panel,  Since  the  vortices  in  the  eight 
neighboring  panels  arc  excluded  from  the  Laurent  series,  simple  geometry  shows 
that  the  convergence  ratio  is  at  least  3/^/2.  Tlierefore.  truncating  the  Laurent  scries 
at  22  terms,  each  vortex  would  be  summed  to  a  relative  erro  6-  10“®;  about  the 
machine  accuracy  in  half  precision  on  the  CYBER  205. 

The  last  step  in  the  procedure  in  Fig.  5  is  the  evaluation  of  the  Laurent  scries  of 
the  curient  inothcr  panel.  While  this  Laurent  scries  could  be  found  using  Eq.  (3b), 
It  can  be  found  more  cfTicicntly  from  the  Laurent  scries  of  the  subpancls.  The 
contribution  of  each  of  the  subpancls  to  the  Laurent  scries  of  the  mother  is  given 
by 
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for  the  first  through  the  fourth  subpancls,  respectively.  {The  factor  2  ‘  m  the  above 
expressions  reflects  the  scaling  of  Laurent  series  proportional  to  the  panel  size.)  The 
coefTicicnts  ,  can  be  evaluated  a  priori  after  which  the  evaluation  of  the 
coefficients  vectorizes. 


4.  PiatioKMANcn 

The  numerical  performance  of  the  prc.scnt  algorithm  is  difl'icult  to  analyze  in 
general.  In  the  following,  the  analysis  has  been  simplified  by  assuming  that  the  ;V 
vortices  arc  homogcncmisly  distributed  over  a  square.  In  that  case,  the  domain  will 
be  subdivided  in  panels  containing  the  s.iinc  number  of  vortices,  n.  each.  The 
number  of  un.subdividetl  panels  is  /V/n  and  the  number  of  levels  of  subdivisions 
needed  is  logjf/V/n). 


The  total  computational  time  fo  find  the  velocity  of  the  ;V  vortices  consists  of  a 
number  of  contributions  First  of  all,  the  vortices  mtist  be  gathered  into  panels.  The 
first  subdivision  of  the  total  domain  involves  -V  vortices  which  are  first  examined  on 
.v-position,  then  on  i-position,  and  correspondingly  reordered.  The  time  involved  in 
this  step  will  be  written  as 


The  factor  Vc,  is  the  lime  needed  to  compare  the  v-  and  r-positions  of  a  vortex  with 
those  of  the  center  of  the  panel,  pass  the  vortex  4  times  through  the  vector  function 
QSVC'Ml’RS  (in  the  205  implementation,  otherwise  to  store  the  vortex  twice),  and 
add  one  to  the  number  of  vortices  in  first  the  right  half  of  the  panel  and  then  to 
the  number  of  vortices  in  the  subpancl,  using  Q^SCN  T. 

The  penalty  factor  /,,  v  expresses  the  overhead  performed  which  is  independent 
of  the  number  of  vortices  involved,  such  as  computing  and  storing  the  panel  infor¬ 
mation  for  the  four  panels  and,  on  the  205,  .^tarting  up  the  vector  operations.  For 
a  'arge  number  of  vortices,  the  operations  for  the  individual  vortices  dominate  the 
total  time  and  /f;  v  will  approach  unity.  However,  for  vector  processors  such  as  the 
205,  the  vector  operations  for  the  individual  vortices  arc  performed  with  such  a 
speed  that  ,v  becomes  apprccia.blc  when  the  number  of  vortices  becomes  less  than 
a  few  hundred.  (I'or  the  simple  vector  operations  in  half  precision  on  the 
FSU/D()E205,  the  start-up  overhead  becomes  equivalent  to  the  time  of  execution 
when  the  number  of  vortices  is  200).  The  subscript  .V  in  the  penalty  factor  /(;  ,v 
expresses  the  representative  number  of  elements  or  vector  length. 

In  the  next  level  of  subdivision,  four  panels  with  each  \N  vortices  are  subdivided, 
requiring  a  computational  time 

/f;  ,V/4  -4  =  .VC,;/(;  V  .: 

Since  there  arc  log,,(A'//M  levels  of  subdivisions  and  the  penally  factor  increases  with 
decreasing  vector  length,  the  total  lime  for  finding  the  panels  may  conservatively  be 
written  as: 

0,  = 'V(-, log,  (9a) 

The  logarithmic  factor  may  be  bounded  by  the  maximum  number  of  subdivisions 
allowed  by  the  machine  ,iccuracy  [8],  but  such  a  bound  depends  on  the  particular 
coding  techniques  ;ind  nuichine  accuracy  available  and  will  be  avoided  here. 

In  the  present  algorithm,  the  original  sum  in  Eq.  (2)  is  used  to  evaluate  the 
\elocity  induced  by  the  n  vortices  in  each  unsubdividcd  panel  upon  its 
neighborhood  of  nine  panels.  If  u^-  is  the  time  needed  to  evaluate  a  single  term  in 
Eq.  (2),  the  total  time  can  be  written,  conservatively,  as 
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neglecting  panels  that  may  fall  outside  the  domain,  or 
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after  correction.  Since  n  will  typically  be  si/'ably  smaller  than  /V,  Eq.  (9b)  will  be 
used. 

For  each  of  the  unsiibdividcd  panels.  A'  cocfricicnts  of  the  Laurent  scries  (3b) 
must  be  determined,  requiring  a  time 


N 

t^=nki\lc,„  (9c) 


These  Laurent  scries  are  ne.xt  used  to  determine  the  velocity  induced  on  the  4  -9/1 
vortices  in  the  neighborhood  of  the  larger  panel  containing  the  unsubdivided  panel, 
e.xcluding  the  9/i  vortices  in  the  neighborhood  of  the  unsubdivided  panel  itself.  The 
time  needed  for  N;n  unsubdivided  panels  is 
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n 

Similarly,  the  Laurent  scries  for  the  /V/4/i  larger  panels  arc  used  to  evaluate  the 
velocity  induced  upon  27  4»  vortices,  neglecting  edge  effects.  With  logj(/'//4)  levels 
of  subdivision,  the  time  can  be  written  conservatively  as: 


N 

!,  =  K-lln  v  J,^,,  -  log4 
n 


(9d) 


The  procedure  of  Grccngard  and  Rokhhn  [7]  avoids  the  logarithmic  factor, 
since  the  number  of  cr'cnicicnts  in  the  Taylor  series  tiocs  not  increase  with  panel 
size.  However,  the  logarithmic  factor  in  (9a)  would  remain. 

Time  is  further  needed  to  combine  the  Laurent  scries  of  the  unsubdivided  panels 
into  these  of  the  larger  panels.  In  the  205  implementation,  each  coefficient  C^.  of  the 
larger  panel  was  written  as  an  inner  product  between  the  vector  of  4A^  coefficients 
of  the  subpanels  and  a  corresponding  vector  of  coefficients  m*,  Eqs.  (8a)  through 
(8g).  The  time  needed  is: 


h 


4  K- 


3 


—  I’/  //.JA  A'- 


(9e) 


The  contribution  most  difficult  to  estimate  is  the  overhead  involved  in  addressing 
the  panel  structure.  For  each  panel,  the  neighborhood  needs  to  be  established,  as 
well  as  the  neighborhood  of  the  next  larger  panel.  Most  of  the  operations  involved 
will  roughly  be  proportional  to  the  number  of  panels:  N/n  unsubdivided  ones,  N/4n 
next  larger  ones,  and  so  on,  a  total  of  less  than  4A'/3/(  panels.  However,  the  binary 


search  to  find  the  nine  neighboring  panels  requires  operations  proportional  to 
log,(4/V/3n).  Tlie  time  for  overhead  will  therefore  be  written  as 
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(9f) 


where  .r„  and  s,  are  representative  computational  times  for  each  panel  and  for  each 
binary  search,  respectively,  neglecting  logj  4;3.  The  symbol  .v  was  used  here  instead 
of  V  in  order  to  indicate  that  the  operations  involved  are  largely  scalar. 

Using  these  various  contributions  to  the  computational  time,  the  decision  when 
to  stop  subdivision  of  the  panels  can  be  addressed.  Collecting  all  contributions,  the 
total  time  needed  to  find  the  velocity  becomes; 

I  =  1  +  27vJ)_  „K+l  -]  N  log,  - 

\  i  II J  II 


+  9i-. /.,,/!  +  Ucfc..,K+  -  vj,  —  + 

\  3  n 


K-  4  .V, 


3  n 


.V. 


:io) 


In  estimating  the  relative  importance  of  the  terms,  it  will  be  assumed  that  the 
nu’  er  of  vortices  N  is  large.  Indeed,  the  number  of  \'orticcs  must  be  sizably  larger 
than  the  typical  mimber  of  terms  in  the  Laurent  senes  in  order  for  the  algorithm 
to  be  useful. 

Under  the  limiting  process  where  both  N  and  the  number  of  vortices  per  panel 
r.  tend  to  infinity,  corresponding  to  relatively  few  large  panels,  the  dominant  term 
m  Eq.  (10)  is  the  time,  Eq.  (9b),  for  the  original  sum,  as  could  be  expected,  oince 
this  term  is  proportional  to  /i,  decreasing  the  number  of  vortices  per  panel  leads  to 
corresponding  reductions  in  computation  time. 

How'cvcr,  when  decreasing  the  number  of  panels,  adverse  affects  must  eventually 
occur  The  penalty  factor  / v  „  increases  when  the  value  of  //  decreases,  since  the 
start  up  time  increases  in  relative  importance.  On  a  tsvo  pipe  205,  the  vector  start 
up  becomes  dominating  when  the  number  of  vortices  becomes  less  than  200, 
limiting  further  reductions  in  the  time  needed  for  the  original  sum. 

On  the  other  hand,  the  time  needed  for  other  operations  increases  while  n 
decreases  For  example,  the  time  for  doing  the  Laurent  series  (9d)  increases  when 
n  decreases  below  a  certain  limit,  since  the  penalty  factor  increases.  This  term 
contains  the  relatively  large  numerical  factor  21K,  so  that  appreciable  increases  in 
the  penalty  factor  tend  to  be  important.  In  addition,  the  scalar  times,  which  can  be 
relatively  large  on  a  205,  arc  inversely  proportional  to  n. 

It  may  be  concluded  that  for  su(Ticicnlly  many  vortices,  the  computational  time 
first  decreases  with  the  number  of  vortices  per  panel  and  then  increases.  As  .i  result, 
a  number  of  vortices  per  panel  c.xists  for  which  the  present  algoritlmi  performs 
optimally 

For  that  reason,  m  generating  the  panels,  the  pieseni  algorithm  decides  whether 


to  subdivide  panels  further  based  on  the  number  of  vortices  in  the  panel.  Further 
subdivisions  are  only  made  when  the  number  of  vortices  is  greater  than  some 
minimum  value  chosen  a  priori. 

Table  I  provides  examples  of  the  influence  of  the  value  of  n  on  the  computational 
time.  In  this  case,  the  vortices  were  approximately  homogeneously  distributed  over 
the  interior  of  a  circle,  grouped  in  rings.  It  appears  from  Table  I  that  the  minimum 
number  of  vortices  to  subdivide  a  panel  on  a  205  should  be  roughly  200. 
Fortunately,  the  precise  value  used  appears  to  have  rclalr  cly  little  influence  on  the 
results. 

In  addition  to  the  computational  time,  the  numerical  errors  in  the  algorithm  arc 
important.  For  p  vortices  of  strength  f'  located  on  a  ring  of  radius  R,  the  velocity 
induced  is 


1 


/  ♦ 

/ 


=  7-> 


pr  zf 

InZ,  Z';-R'’ 


TABLE  I 


Compu..i(ion.il  Time  .iml  Numerical  Errors  Tor  Vorliccs  Homo('ciiccnislv  Dislrihulcd 
wiihm  a  Circle  Using  2.1  Term  Asymptotic  lixpaiisions 


Number 

of 

vornccs 

1000 

:(HH) 

4000 

XlH)0 

1  6,(X)0 

12.IX)0 

64.000 

Time  for  summation,  CPI' 

.sCC(>IHis 

Original 

0.10 

0.10 

I.1S 

5  19 

2t  15 

X6  12 

156  IS 

IPO 

0,11 

0  >5 

0.71 

1  95 

' 

0  u 

17.17 

200 

O.ll 

0  27 

067 

161 

s  12 

17  10 

4iV1 

0.10 

027 

OS2 

1  61 

■1 

S  12 

5.1 

Ratio  of  improvement 

21X1 

0.9 

1  1 

20 

yy 

6  0 

10  6 

M.iMmum  ciior  in  Ihc  vcIt* 

city,  percent 

Original 

0.005 

0019 

0.020 

0040 

0  oso 

0  1 60 

n 

nx) 

0004 

otxn 

0(X19 

0  013 

0016 

0  022 

11024 

21XJ 

0.004 

0  IX)f> 

0(X)9 

0O11 

0016 

0  024 

41)0 

0,005 

O'XIO 

0012 

0041 

0  021 

0  021 

0  010 

Mean  square  error  in  the  'clocili,  pcicent 

Original 

0.00.1 

01X16 

0012 

0021 

ti  046 

0  1S7 

200 

0.002 

01X14 

0  005 

0iX!7 

0.009 

0  012 

0  011 

Average  error  i 

n  the  vcloctt 

y.  percent 

Original 

0.001 

0010 

0  020 

0041 

0  0S2 

0  1 64 

200 

0,002 

DimU 

0  005 

0(X)6 

I'OOS 

01)10 

0  012 

-ri 


or 


=  V  -  I 


{p-\)r 

4nZ, 


if  the  vortex  /  is  located  on  the  same  ring.  By  summation  over  all  rings  the 
analytical  solution  can  be  found  and  compared  to  the  obtained  results. 

Table  I  list  the  maximum  deviation  in  cither  velocity  component  from  the 
analytical  solution,  expressed  in  a  percentage  of  the  velocity  at  the  perimeter  of  the 
circle.  The  present  algorithm  shows  considerably  better  accuracy  than  the  original 
sum,  which  may  be  due  to  the  summation  of  the  terms  in  groups.  In  the  original 
algorithm,  the  individual  terms  were  added  to  an  increasingly  large  total,  leading 
to  a  loss  of  signincant  digits. 

i-'or  random  walk  computations,  the  mean  square  or  average  errors  may  be  more 
relevant  than  the  maximum  error,  since  only  .iveraged  quantities  arc  relevant.  Both 
these  errors  show  behavior  similar  to  the  maximum  error. 


TABl.n  II 


A'  Fable  I.  I'.il  Using  13  Term  Expansions 


Number 

of 

vorliecs 

ItXX) 

;ixK) 

4(XXf 

sooo 

1 6,000 

32,tX10 

64,000 

Time  for  suinmalion,  Cl’tJ  seconds 

Original 

0,10 

0  36 

I.3K 

5.39 

21  36 

R6,35 

356,72 

100 

0.09 

0  26 

0.57 

1  40 

2  79 

6  61 

13,03 

;oo 

0,09 

It  :5 

0  55 

UI7 

2  79 

691 

1 3  08 

400 

0.10 

0  25 

0  79 

1  43 

4  27 

6,92 

20,40 

Ratio  of  improv.emenl 

100 

1.1 

1  4 

24 

3  9 

7  7 

13.1 

27  4 

M.iximuin  error  in  the  velocity 

percciU 

Original 

0.005 

0010 

0.020 

0  040 

0.080 

0.160 

0.319 

100 

0.003 

0  005 

0.007 

0  009 

n.oi  1 

0.013 

0  015 

:oo 

0.003 

0006 

0.007 

0011 

0.0 1 1 

0.016 

0.015 

400 

0.005 

0  006 

0.01 1 

0012 

0.019 

0016 

0  025 

Mean  square  error  in  tlic  vclociiy.  percent 

Original 

0003 

0  006 

0.012 

0023 

0.046 

0,093 

0  187 

100 

0.002 

0  003 

0004 

0  005 

0.006 

0007 

0  009 

Average  crrc'r 

in  the  velocity,  j 

vreent 

Original 

0.003 

u(X)5 

0  010 

0  020 

0  04 1 

00S2 

0  164 

100 

0.002 

tUX)2 

0003 

otXM 

0  005 

0,006 

0  008 

> 


The  results  in  Table  I  were  obtained  by  expanding  all  Laurent  series  to  machine 
precision,  23  terms.  Yet  in  view  of  the  final  errors  in  the  results,  there  appears  little 
justification  in  demanding  such  accuracy,  unless  special  provisions  arc  made  to 
avoid  accumulation  of  the  round-off  errors.  Results  for  13  term  Laurent  scries  arc 
presented  in  Table  II.  Remarkably,  the  resulting  errors  prove  somewhat  lower  than 
those  in  the  23  term  expansion.  A  good  explanation  of  this  effect  cannot  be  given; 
however,  the  maximum  possible  error  in  truncating  the  Laurent  series  is  only 
O.CK)6%,  small  compared  to  the  final  errors.  On  the  other  hand,  the  final  terms  in 
the  Laurent  series  correspond  to  the  fastest  l-'ouricr  components:  for  that  reason 
truncating  the  series  may  have  some  averaging  cfTcct  on  the  round-off  errors. 

For  the  case  of  Tables  I  and  II.  the  vorticity  occupied  most  of  the  domain  under 
consideration.  A  somewhat  different  :asc  arises  when  the  vortices  are  evenly  spaced 
along  the  perimeter  of  a  circle.  .Since  the  vorticity  is  now  sparsely  distributed,  the 
present  procedure  will  generate  a  considerable  number  of  empty  panels,  and  it 


TAIILR  III 


Compulatioii.Tl  l  inic  ami  Numerical  Errors  for  Vorliccs  Homogeneously 
Dislribulecl  on  a  Circle  Usinp  13  Term  Asymplolic  Expansions 


Number 

of 

Noriiccs 

1000 

:iKX) 

4u00 

SOOO 

16,000 

32,000 

64,000 

Time  for  summation,  Cl’ti 

seconds 

Original 

0,10 

It  .>6 

1  38 

541 

21  38 

86  32 

356,75 

50 

008 

0  17 

0  36 

0,78 

1  65 

3  53 

7.30 

100 

0,06 

0  14 

0.3t 

0  65 

1  42 

3  04 

6.28 

:iX) 

0,06 

0  17 

0  34 

0  75 

151 

3  36 

7.19 

Raiio  of  improvement 

100 

1,7 

2 

4  5 

83 

15.1 

28  4 

55.8 

Maximum  cror  in  the  velocity,  percent 

Ongjnal 

0.02-) 

OU6I 

0  128 

0  271 

0  547 

1  100 

— 

50 

0012 

0  023 

0IW6 

0093 

0.180 

0  360 

0.702 

100 

0  012 

0  024 

0  046 

0  094 

0,181 

0,361 

0.702 

:oo 

0,013 

0  026 

0049 

0  098 

0  1 84 

0  363 

0.705 

Me.in  square  error  ,n  the  vclocilx,  percent 

Orijin.il 

0.005 

0  010 

0018 

0  041 

0.079 

0  143 

— 

100 

O.OOA 

0  008 

0  012 

0  031 

0059 

0096 

0.149 

Average  crroi 

in  llic  vclociiy.  pcrccnl 

Orijin.il 

0  001 

0  002 

0003 

0  006 

0012 

0  024 

— 

100 

otxw 

0001 

0001 

0  001 

t)n02 

0  002 

0  002 

/  ^7 
/  ■ 


niighl  sccin  that  this  would  adversely  affect  coinpiilalional  lime.  However,  ihe  ilata 
in  Table  III  shows  that  performance  improves.  This  appears  to  agree  with  the 
observations  of  Carrier.  Greengard.  and  Rokhiin  [8]. 

The  simplined  ease  of  vortices  arranged  on  a  single  horizontal  line  can  shed  some 
light  on  this  increase  in  cfTiciency.  Repeating  the  previous  analysis  with  suitable 
modifications,  the  total  time  becomes: 


'  —  2a,,  / ,,  „  +  Cl',  f,  „ K  +  8  — 

/I 


+  +  >-'cfc.n^  + i,.iK - 1-4—  N. 

n  n 


Comparison  with  Eq.  1 10)  does  show  that  the  time  for  scalar  panel  overhead  has 
increased.  On  the  other  hand,  the  time  needed  for  the  direct  summation  has 
improved,  since  the  empty  panels  decrease  the  number  of  vortices  in  the 
neighborhood  of  unsubdivided  panels  from  9n  to  3/i.  In  addition,  the  time  for  the 
Laurent  series  has  decreased,  rmcc  6/i  rather  than  27n  vortices  need  to  be 

TAnLFZ  IV 

('ompiilaiioniil  Time  and  Numerical  Errors  for  Voriiec,'; 
on  the  Perimcicr  of  a  Circle  for  a  MieruVAX  II 


Number 

of 

voriices 

40(1 

soo 

1600 

3200 

6400 

12,800 

25,600 

Time  for  oi 

mmalion,  CPU  .seconds 

Origin.il 

9,1 

36,3 

151  9 

614,2 

2470  X 

9867,3 

39469,2 

30 

5,0 

109 

24,5 

56.3 

1270 

272  3 

604.1 

40 

4,7 

106 

25  1 

554 

128.5 

286  3 

597.8 

50 

4,7 

106 

25,0 

55.5 

1.33  5 

292.5 

616.3 

60 

4,7 

10,6 

26,4 

55.7 

138.9 

303.3 

635  6 

Mean  square  error  in  llic 

veloeny,  perccnl 

Original 

OOOOtKI 

lUXXXX) 

0  01X301 

0  (XX304 

0,00007 

0.00019 

0.00036 

40 

Q.OOOOl 

(lOOOOl 

0(XX)0t 

0  00002 

0,tX)004 

0  00010 

OOOOIS 

.Addiiional 

.irr.ny  storage 

used.  4  b)ie  rsords 

40 

496 

3:S9 

.3907 

5040 

7281 

1 1 374 

20827 

Pcrccnlagc 

CPU  lime  for 

>fcps 

62 

50 

4S 

43 

42 

41 

34 

U, 

30 

14 

S 

9 

7 

7 

0: 

0 

26 

34 

40 

4  3 

50 

determined  from  the  Laurent  scries  of  the  unsubdivided  panels;  the  remaining 
contributions  are  found  using  the  combined  Laurent  series. 

Next,  while  scalar  operations  are  relatively  slow  on  the  205,  the  corresponding 
total  times  are  proportional  to  l/n,  which  is  of  order  10"^.  On  the  other  hand,  the 
times  for  the  direct  sum  and  the  Laurent  scries  are  proportional  to  9n  and  UK, 
respectively,  which  are  of  order  10^.  For  that  reason,  the  scalar  times  need  not 
dominate  even  for  quite  slow  scalar  processing. 

While  the  present  algorithm  still  leads  to  reduced  errors,  the  differences  are  not 
so  pronounced  as  in  Tables  I  and  II.  The  reason  may  be  that  in  this  case  the 
vortices  in  the  immediate  vicinity  of  each  other  dominate  the  errors  (the  sum 
approximates  a  singular  integral).  Those  vortices  arc  still  summed  in  the  same  way 
as  in  the  original  algorithm. 

Table  IV  shows  computational  times  for  a  scalar  version  of  the  algorithm.  As 
may  be  expected,  the  scalar  version  can  address  smaller  groups  of  vortices  more 
efficiently  than  the  205  version,  resulting  in  some  additional  savings. 

The  total  time  CPU  time  used  can  be  divided  into  contributions  to  to  perform  the 
original  sum  (2),  ,  to  perform  the  Laurent  series  of  panels  which  are  not  further 

subdivided,  and  i , ,  for  the  Laurciu  scries  of  panels  which  arc  Table  IV  shows  those 
contributions  for  the  ctise  that  n  --40.  Recasting  of  Laurent  scries  into  Taylor  series 
as  used  by  Grccngard  and  Rokhlin  [7]  could  be  used  to  reduce  the  required 
time 

Comparison  of  Table  IV  with  the  data  of  Carrier.  Grccngard,  and  Rokhlin  [8] 
does  suggest  a  significant  increase  in  storage  due  to  such  a  recasting. 


5.  CoNCl.fDING  RhMARKS 

The  present  algorithm  is  concerned  with  fast  solution  of  the  2-dimcnsiotial 
F’oisson  equation  under  pointwise  forcing.  Since  the  actual  application  is  not  the 
true  subject  of  this  paper,  or.'  a  concise  description  of  the  one  considered  here  will 
be  given.  Lagraiigian  flow  computations  using  a  random  walk  simulation  of 
diffusion  effects  similar  to  [9].  A  relatively  simple  removal  of  the  singular  behavior, 
Eq.  (4),  was  used  in  computed  examples  such  as  Fig.  I.  In  most  computations,  the 
chosen  vortex  diameter  was  0  675  times  the  random  step  si/.c  n  =  riic 

normal  boundary  condition  was  satisfied  by  means  of  mirri'r  \ortices  within  the 
cylinder.  To  satisfy  the  l.ingcntial  boundary  condition,  after  each  predictor-correc¬ 
tor  step  all  vortices  within  a  distance  1.27  rr  from  the  wall,  a  thin  sub-lavcr  of  the 
boundary  layer,  were  removed.  Next  the  slip  velocity  at  the  wall  was  evaluated  and 
integrated  to  find  the  amount  of  circulation  needed  to  satisfy  the  no-slip  condition 
This  circulation  was  subsequently  assigncil  to  a  ring  of  vortices  a  distance  0.675  n 
away  from  the  wall  .iiui  spaced  1.27  rr  apart.  I'hc  procedure  leads  to  the  same  flux 
of  vortices  through  the  cutoff  at  I  27  a  as  a  homogeneous  distribution  of  vortices 
within  the  cutoff.  In  order  to  reduce  the  random  fluctuations  introduced  by  strong 
xortices,  the  number  of  \orticcs  placed  at  each  location  along  the  wall  was  chosen 


10  give  an  approximately  uniform  vortex  strength.  At  least  one  vortex  was  placed 
at  each  location  if  the  local  circulation  was  non-zero.  Some  experiments  varying  the 
given  numerical  values,  or  using  an  exponential  vortex  core  rather  than  Eq.  (4), 
were  performed,  but  results  were  ambiguous  due  to  the  random  noise.  The  random 
step  sizes  were  taken  from  a  data  base  of  8000  random  numbers,  starting  from  a 
randomly  chosen  position. 

The  purpose  off  this  paper  was  to  sliow  how  the  computational  time  can  be 
greatly  reduced  using  Laurent  scries,  allowing  a  much  larger  number  of  vortices  to 
be  included.  The  use  of  Laurent  scries  or  replacement  elements  to  save  computa¬ 
tional  time  is  not  a  new  notion  [6];  however,  the  present  method  renders  the 
application  effective  by  gathering  the  point  forces  into  an  adaptive,  ordered  panel 
structure.  The  contribution  of  the  present  paper  is  therefore  primarily  a  program¬ 
ming  technique  which  allows  an  easily  addressable  adaptive  description  of  irregular 
distributions  of  points.  Moreover,  it  is  quite  suited  for  vector  processing  and 
requires  little  storage.  It  seems  simpler  and  possibly  more  vcctorizable  than  the 
procedure  of  Carrier,  Grccngard,  and  Rokhlin  [8]. 

The  evidence  of  the  Tables  1  through  IV  shows  that  the  present  algorithm  is 
"fast”  in  the  sense  that  the  computational  time  roughly  doubles  when  the  number 
of  vortices  doubles.  For  the  original  sum  in  Eq.  (2),  the  computational  time 
becomes  larger  by  a  factor  four  instead.  For  that  reason  the  savings  in  computa¬ 
tional  time  increase  with  the  number  of  vortices. 

In  fact,  the  time  estimates  in  Eqs.  (10)  and  (12)  show  the  computational  time  to 
be  proportional  to  /VIogA',  similar  to  the  fast  Fourier  transform  solutions  of  the 
Poisson  equation  such  as  Hockney's  FACR  algorithm,  which  needs  Alog2(log2  A) 
operations.  Mowever,  a  closer  study  of  Eqs.  (10)  and  (12)  shows  that  for  typical 
valuesK~20  and  /i~100.  the  coefficient  of  the  /VIogA  term  in  the  present 
algorithm  will  be  numerically  quite  large. 

For  that  reason,  one  of  the  motivations  mentioned  m  the  Introduction  should 
still  be  present  in  order  to  adopt  an  algorithm  such  as  the  present  one. 

.■Ml  interesting  question  is  whether  the  present  algorithm  is  applicable  to 
3-dimcnsional  Poisson  problems.  This  would  make  it  possible  to  address  such 
problems  as  the  motion  of  stars  in  galaxies  and  3-dimcnsional  flows  with  sparse 
vortex  geometry.  Most  of  the  procedures  in  Sections  2  and  3  carry  through 
immediately  by  the  simple  step  of  including  the  digits  of  the  third  coordinate  in  the 
panel  numbers.  However,  the  straightforward  generalization  of  the  Laurent  series  to 
spherical  harmonics  as  applied  by  Greengard  and  Rokhlin  [II]  has  the  disadvan¬ 
tage  that  the  number  of  terms  added  for  each  order  of  accuracy  increases  A  proce¬ 
dure  based  on  fast  Fourier  transforms  proposed  by  Greengard  and  Rokhlin  [12] 
can  significantly  reduce  the  effort. 

The  present  procedure  of  generating  and  addressing  a  complex  panel  structure 
docs  not  need  to  be  restricted  to  solution  of  the  Poisson  equation,  but  could  be 
used  for  other  problems  involving  groups  of  points  in  which  (he  interaction  between 
elements  of  different  groups  can  be  simplined  when  the  distance  between  the  groups 
is  sufficient. 
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APPENDIX  III 

Summary  of  the  Computational  Scheme. 


J 


Computations 


Random-walk  vortex  simulations  of  the  full  Navier-Stokes  equations  were  performed 
as  comparison.  In  the  computations,  the  flow  field  was  represented  by  discrete  vortex  blobs 
of  the  form: 

_ Vj  —  Vi _ 

{xj  -  Xi)^  +  {yj  -  yi)2  -I-  d? 

_  _  Xj  —  Xj _ 

i  ~  +  (yj  ~ 

Ti 

A  recently  developed  fast  procedure  allowed  the  computations  to  be  performed  with  sizably 
more  vortex  blobs  than  previously  possible  [see  J.  Comp.  Phys.  paper].  The  vortices  were 
advanced  in  time  using  a  two-step  Runge-Kutta  scheme.  To  simulate  diffusion  effect, 
each  time  step  the  vortex  motion  was  augmented  with  a  random  component  of  average 
magnitude 

The  normal  wall  boundary  condition  was  satisfied  by  mirror  vortices,  after  a  mapping 
of  the  airfoil  onto  a  circle.  The  mapping  used  was  a  generalized  Von  Mises  transform  which 
correctly  reproduces  the  kinks  in  the  contour  at  the  trailing  edge: 


dZ 


~ih 


The  constants  (7,  and  7a,  are  determined  from:  the  Idnks  in  the  NACA  0012  airfoil 
contour  caused  by  the  small  but  finite  thickness  at  the  trailing  edge,  the  regularity  of  the 
mapping  at  infinity,  and  finally  from  least  square  minimization  of  the  errors  in  airfoil  shape 
elsewhere.  The  Von  Mises  type  procedure  was  preferred  above  a  Fast  Fourier  transform, 
since  the  transform  is  relatively  inaccurate  due  to  the  singularities  in  contour.  In  addition, 
the  transform  would  be  quite  inefficient  during  the  actual  flow  computation. 

The  no-slip  boundary  condition  was  satisfied  by  addition  of  vortices  at  the  wall  during 
each  time-step:  First  all  vortices  within  a  distance  of  1.27\/2r/Af  were  removed.  Then  a 
ring  of  new  vortices  was  added  at  a  distance  0.675  \/2i>' At  to  correct  the  wall  velocity  to  zero. 
(The  distance  for  sidding  vortices  equals  the  diffusion  distance  of  the  vorticity  generated 
by  the  wall  during  the  time-step  for  the  true  Navier-Stokes  equations;  the  removal  distance 


was  chosen  based  on  a  statistical  study  requiring  that  the  scheme  handles  locally  uniform 
vortirity  distributions  accurately,  not  unlike  discretization  technifines  in  fitiite  difTerence 
procedures).  The  vortex  diameter  dj  was  rather  arbitrarily  chosen  to  be  0.675\/2i/At; 
testing  showed  that  results  depended  little  on  the  actual  value  used. 

The  CYBER  205  results  were  post-processed  on  the  departmental  MicroVAX  II,  using 
a  fast  Fourier  transform  to  find  the  streamlines.  The  vorticity  was  represented  in  bit¬ 
mapped  graphics  as  half  tones. 
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Abstract 


Although  unsteady,  high-Reynolds-number,  laminar  boundary  layers  have  convention¬ 
ally  been  studied  in  terms  of  Eulerian  coordinates,  a  Lagremgian  approach  may  have  sig¬ 
nificant  analytical  and  computational  advantages.  In  Lagrangian  coordinates  the  classical 
boundary-layer  equations  decouple  into  a  momentum  equation  for  the  motion  paredlel  to 
the  boundary,  and  a  hyperbolic  continuity  equation  (essentially  a  conserved  Jacobian)  for 
the  motion  normal  to  the  boundary.  The  momentum  equations,  plus  the  energy  equation  if 
the  flow  is  compressible,  can  be  solved  independently  of  the  continuity  equation.  Unsteady 
separation  occurs  when  the  continuity  equation  becomes  singular  as  a  result  of  touching 
characteristics,  the  condition  for  which  can  be  expressed  in  terms  of  the  solution  of  the  mo¬ 
mentum  equations.  The  solutions  to  the  momentum  and  energy  equations  remain  regular. 
Asymptotic  structures  for  a  nmnber  of  tmsteady  three-dimensional  separating  flows  follow 
and  depend  on  the  symmetry  properties  oi  the  flow  (e.g.  line  symmetry,  axial  symmetry). 
In  the  absence  of  any  symmetry,  the  singularity  structure  just  prior  to  separation  is  found 
to  be  qusisi  two-dimensional  with  a  displacement  thickness  in  the  form  of  a  crescent  shaped 
ridge.  Physically  the  singularities  can  be  understood  in  terms  of  the  behavior  of  a  fluid 
element  inside  the  boundary  layer  which  contracts  in  a  direction  parallel  to  the  boundary 
and  expands  normal  to  it,  thus  forcing  the  fluid  above  it  to  be  ejected  from  the  boundary 
layer. 


1.  Introduction 

A  major  feature  of  unsteady  large-Reynolds-number  flow  past  a  rigid  body  is  the  shed¬ 
ding  of  vortices  from  the  surface  of  the  body.  Such  vortices  alter  the  forces  exerted  on  the 
body  dramatically  (McCroskey  &  Pucci  1982).  A  more  complete  theoretical  understanding 
of  vortex  shedding  would  be  advantageous  in  the  design  of  air,  land  and  water  transport. 
Theoretical  models  of  vortex  shedding  also  have  application,  inter  alia,  in  the  description 
of  air  flow  over  hills  and  water  waves,  water  flow  over  sand  ripples,  and  blood  flow  through 
curved  and  constricted  arteries  and  veins. 

A  classical  example  of  vortex  shedding  develops  when  a  circiJar  cylinder  is  set  into 
motion  in  the  direction  normal  to  its  axis.  This  example  was  first  studied  by  Prandtl  (1904), 
and  the  process  by  which  an  initially  attached  boundary  layer  develops  into  a  separated  flow 
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with  detached  free  shear  layers  has  been  clearly  illustrated  by  the  experiments  of  Nagata, 
Minami  &  Murata  (1979),  and  Bouard  &  Coutanceau  (1980).  The  term  ‘separation’  will 
in  this  paper  be  used  to  refer  to  the  ‘breakaway’  of  a  thin  layer  of  vorticity  from  the 
surface  of  a  body.  This  definition  of  separation  is  close  to  that  of  both  Prandtl  (1904) 
and  Sears  &  Telionis  (1975).  In  particular,  Sears  &  Telionis  speak  only  of  separation  when 
the  penetration  of  the  boimdary-layer  vorticity  away  from  the  wall  becomes  too  large  to 
be  described  on  the  usual  0{Re~*)  boundary -layer  scale  {Rt  is  the  Reynolds  number  of 
the  flow,  and  is  assumed  Iwge).  Therefore,  once  separation  has  developed  the  classic^J 
attached  flow  solution  will,  in  general,  no  longer  be  valid. 

The  first  theoretical  advance  in  understanding  the  unsteady  cylinder  flow  at  high 
Reynolds  numbers  was  made  by  Blasius  (1908).  He  explained  the  occurrence  of  flow 
reversal  inside  the  attached  unsteady  boundary  layer  which  is  set  up  immediately  the 
cylinder  starts  to  move.  In  the  case  of  steady  flow  past  a  rigid  surface,  flow  reversal  is 
often  accompanied  by  separation.  However,  Moore  (1958),  Rott  (1956)  and  Sears  (1956) 
all  realized  that  zero  wall  shear  is  not  necessarily  related  to  separation  in  unsteady  flow. 
Sears  &  Telionis  (1975)  noted  subsequently  that  their  definition  of  separation  is  consistent 
with  the  termination  of  the  boundary -layer  solution  in  a  angularity.  Such  a  angularity 
will  be  referred  to  as  the  separation  singularity,  and  the  time  at  which  it  develops  as  the 
separation  time. 

A  considerable  number  of  numerical  computations  have  attempted  to  verify  the  ex¬ 
istence  of  a  singularity  in  the  boundary-layer  solution  for  the  circular  cylinder  problem. 
The  first  convincing  evidence  that  a  singularity  forms  within  a  finite  time  was  given  by 
Van  Dommelen  &  Shen  (1977,  1980a,  1982).  In  a  Lagreuigian  computation,  with  fluid  par¬ 
ticles  as  independent  coordinates,  they  found  that  a  separation  singularity  develops  after 
the  cylinder  has  moved  approximately  |  of  a  diameter.  The  existence  of  this  singularity 
has  been  confirmed  by  the  finite  difference  numerical  calculations  of  Ingham  (1984)  and 
Cebeci  (1982)  (however  see  Cebeci,  1986),  and  the  computer  extended  series  solution  of 
Cowley  (1983).  These  calculations  were  all  based  on  Eulerian  formulations.  A  similar 
two-dimensional  separation  singularity  has  been  observed  using  Lagrangian  procedures  on 
an  impulsively  started  ellipse  at  several  eingles  of  attack  (Van  Dommelen,  Wu,  Chen  & 
Shen,  unpublished  results),  on  airfoils  (Wu  1988),  in  turbulence  production  (Walker  1988), 
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on  an  impulsively  started  sphere  (Vaui  Dommelen  1987),  and  using  Eulerian  schemes  in 
leading  edge  stall  (Cebeci,  Khattab  &  Schimke,  1983)  and  about  a  rotating  cylinder  (Ece, 
Walker  &  Smith,  1984). 

Ekcluding  vortex  methods,  flows  with  free  surfaces,  and  some  more  specialized  com¬ 
pressible  flow  computations,  Lagrangian  coordinates  have  not  been  as  widely  used  as  their 
Eulerian  counterparts  in  fluid  mechanics,  especisJly  for  boundary-layer  flows.  Yet  for  some 
flows,  such  as  unsteady  flows  in  which  advection  dominates  diffusion,  Lagrangian  coordi¬ 
nates  seem  more  appropriate  (e.g.  see  the  in  viscid  calculations  of  Stern  &  Paldor  (1983), 
Russell  &  Landahl  (1984)  and  Stuart  (1987)).  As  far  as  unsteady  separation  is  concerned, 
the  advantage  of  a  Lagrangian  approach  stems  from  the  fact  that  in  these  coordinates 
the  classical  boundary-layer  equations  decouple  into  a  momentum  equation  for  the  motion 
parallel  to  the  boundary,  and  a  continuity  equation  for  the  motion  normal  to  the  boundary 
(Shen  1978).  The  solution  of  the  former  equation  can  be  found  independently  of  the  latter. 
Moreover,  while  the  time  that  the  separation  singtilarity  develops  can  be  identified  from 
the  solution  to  the  momentum  equation,  only  the  solution  to  the  continuity  equation  is 
singular  (see  section  2). 

An  important  consequence  of  the  Lagrangian  approach  is  that  simple  descriptions  can 
be  found  to  a  wide  variety  of  separations  in  one-,  two-  and  three-dimensional  imsteady 
flows.  In  this  paper  we  consider  unsteady  flows  in  general,  then  in  part  2,  (Van  Dommelen 
1989),  the  separation  process  that  occurs  at  the  equatorial  plane  of  a  sphere  which  is  set 
into  a  spinning  motion  is  examined  in  detail. 

In  the  next  section  we  develop  the  simple  analytic  machinery  needed  to  find  self- 
consistent  three-dimensional  separation  structures  for  both  compressible  and  incompress¬ 
ible  fluids.  Some  of  the  properties  of  the  Lagrangian  version  of  the  boundary-layer  equa¬ 
tions  are  also  discussed.  In  section  3  the  Lagrangian  structure  for  three-dimensional  sep¬ 
aration  is  derived  under  the  assumption  that  the  flow  can  be  completely  general,  then  in 
section  4  the  changes  in  structure  are  discussed  when  various  symmetries  restrict  the  flow 
geometry. 


2.  Lagrangian  Formulation 


The  Lagrangian  description  of  boundary-layer  flow  uses  fluid  particles  (i.e.  infinitesi- 
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mal  masses  of  fluid)  as  the  basis  of  the  coordinate  system.  A  convenient  coordinate  system 
for  the  fluid  particles  is  given  by  the  initial  Eulerian  position  of  the  particles  (see 

Lamb  (1945)  for  example): 


^  at  t  =  0  •  (2.1) 

The  precise  form  of  the  Lagrangian  solution  depends  on  the  particular  reference  time, 
defined  here  as  the  start  of  the  motion,  but  the  physical  solution  is  independent  of  it. 

Following  Rosenhead  (1963)  we  assume  that  the  position  coordinates  x  and  z  describe 
an  orthogonal  coordinate  system  on  the  surface  of  the  body  in  question.  The  lengths  of  the 
line  elements  dx  and  dz  are  taJcen  as  hi  dx  and  h^dz  respectively.  The  coordinate  normal 
to  the  surface  is  denoted  by  j/,  which  is  scaled  with  the  square  root  of  the  reference  shear 
viscosity. 

In  Lagrangian  coordinates,  conservation  of  volume  for  a  compressible  fluid  can  be 
expressed  in  terms  of  a  Jacobian  determinant  as  follows  (e.g.  Hudson  1980): 

pH{x,z)J{x,y,z)  =  poHo  ,  {2.2a) 


where 


J{x,y,z)  = 


^,v 

y,i 

y,v 

^,v 

x,z) 

=  h 

y,< 


(2.26, c) 
{2.2d,  e) 


p{i,r]X,t)  is  the  density  of  the  fluid,  and  a  subscript  comma  denotes  a  Lagrangian  deriva¬ 
tive.  The  velocity  components  of  the  flow  are  related  to  the  fluxions  of  position  by 


u  =  hi{x,z)x  ,  u;  =  63(1,2)2  , 


(2.3a, 6) 


where  a  dot  represents  a  Lagrangian  time  derivative. 

For  compressible  flow,  the  momentum  and  energy  equations  are  (e.g.  Rosenhead 
1963): 

w  I 

p{u  +  {uhlt  -  wh3t)  —  )  =  ~—px  +  DyifiDyu)  +  pgt  ,  (2.3c) 

n  fii 

XL  1 

Piw  +  {whi^  -  uhig)-~)  =  ~—p,  +  Dy{pDyw)  +  pg,  ,  (2.3d) 

M  ^3 
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P^p-^P^'P-^P  =  ti{{DyuY +{DywY)  +  Dy{^DyT)  ,  (2.3e) 

where  /i  is  the  scaled  shear  viscosity,  <t  is  the  Prandtl  number,  and  and  g,  are  the 
components  of  the  acceleration  of  gravity.  The  temperature,  T,  and  internal  energy  e,  are 
assumed  to  be  functions  of  density  and  pressure,  while  the  pressure,  p,  is  a  known  function 
of  z,  z  and  t\  thus 

u  w 

P  =  P.  +  +  ^P.  .  (2.3/) 

For  an  incompressible  flow  p  =  0  and  e  is  taken  to  be  a  function  of  T  and  p. 

Although  the  y-derivative  Dy  is  Eulerian  in  nature,  it  can  be  written  in  the  Lagrangian 
form  (see  also  Shen  1978): 


^  _  p(^,V,Ct)H{x,z)J{x,u,z) 


(2.4) 


FVom  (2.4)  it  follows  that  at  a  fixed  wall  the  Eulerian  Dy  and  Lagrangian  d/dr)  operators 
differ  only  by  the  density  ratio,  which  leads  to  simplifications  in  the  calculation  of  the  wall 
shear. 

Allowing  for  a  moving  boundary,  appropriate  boundary  conditions  to  (2.3)  are: 

(u,ti;,p)  =  («(,(*, z,t),u;j,(x, 2, t),pi(®,z,/))  on  y  =  0  ,  (2.5a) 

{u,w,p) {u^{x,z,t),w^{x,z,t),p^{x,z,t))  as  y oo  ,  (2.56) 

where  Uf,  and  Wf,  specify  the  velocity  of  the  boundary  in  the  x-  and  z-directions  respectively, 
tie  and  We  are  the  corresponding  external  slip-velocities,  pe  is  the  external  flow  density, 
and  the  wall  density  pt,  can  be  given  implicitly  as  the  temperature  at  the  wall.  Ordinarily, 
these  boundary  conditions  translate  immediately  to  the  Lagrangian  domain  by  means  of 
(2.3a,b).  In  the  case  of  suction  or  blowing  through  the  wall,  they  must  be  applied  at  an  t/- 
boundary  moving  through  the  Lagrangian  domain,  however,  the  wall  boundary  conditions 
turn  out  to  be  of  little  importance  for  the  local  analysis  of  this  paper. 

The  principle  advantages  of  Lagrangian  coordinates  derive  from  the  absence  of  both 
the  normal  particle  position  y  and  the  normal  velocity  component  v  from  (2.3)  ariu  (2.4). 
Consequently,  the  particles’  motion,  as  projected  onto  the  surface  of  the  body  (x,z),  can 
be  found  independently  of  the  normal  particle  position  y.  Subsequent  integration  of  the 
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Jacobian  (2.2)  along  lines  of  particles 
particle  position 


at  constant  projected  position  (x,  z)  yields  the  normal 


f 


Po  Hod^ 


pH\Vx  A  Vz\ 


(2.6) 
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where  +  drj^  +  dC^  =  (® Lagrangian  gradient,  and  the 

integral  is  performed  in  the  Lagrangian  (^,n»C;0  coordinate  system  along  the  lines  of 
constant  x,  z,  and  t,  i.e.  lines  which  in  physical  space  are  vertic^J  through  the  boundary 
layer. 

The  central  issue  of  this  paper  can  now  be  stated:  we  hypothesize  that  during  the  evo¬ 
lution  toward  separation,  the  projected  position  (x,z)  can  remain  regular,  and  commonly 
does  remain  regular.  When  true,  such  regularity  strongly  restricts  the  possible  behavior  of 
X  and  z  near  separation,  and  to  characterize  separation  we  need  only  identify  the  nature 
of  solutions  to  the  continuity  equation  (2.2)  or  (2.6)  -  an  equation  which  is  much  simpler 
than  the  momentum  equations.  The  remaining  ambiguity  in  the  behavior  of  x  and  z  is 
resolved  using  arguments  of  symmetry. 

Various  arguments  to  justify  our  hypothesis  can  be  given.  One  of  them  is  self- 
consistency.  If  it  is  assumed  that  x,z,u,w,  and  p  are  non-singular  at  the  separation  time 
t,,  then  the  solution  to  the  Lagrangian  momentum  equations  can  be  expanded  in  powers 
of  (<  —  tf)  to  any  algebraic  order.  In  contrast,  the  usual  Eulerian  asymptotic  expansions 
show  only  that  the  first  few  terms  in  the  expansions  are  self-consistent. 

As  another  argument,  Vain  Dommelen  (1981)  showed  analytically  that  the  inviscid 
incompressible  two-dimensional  equations  have  solutions,  x,  z,  which  are  regular  func¬ 
tions  of  the  Lagrangian  variables,  although  y(4,0  is  singular  (this  analysis  can  be  further 
developed  by  expanding  in  powers  of  a  small  coefficient  of  viscosity).  Yet  this  example  is 
somewhat  artificial;  physically  it  would  require  that  during  the  evolution  of  the  boundary 
layer  the  coefficient  of  viscosity  wets  changed  significantly  by  some  external  means. 

A  more  powerful  eu-gument  is  possibly  the  capability  of  the  analysis  in  this  paper  to 
reproduce  and  extend  several  known  separation  processes  previously  analyzed  in  Eulerian 
coordinates.  However,  the  most  convincing  argument  is  provided  by  actual  numerical  solu¬ 
tions  of  the  Lagrangian  boundary-layer  equations.  For  example.  Van  Dommelen  &  Shen’s 
(1980a)  computation  of  the  boundary  layer  on  aji  impulsively  started  circular  cylinder 
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provided  direct  numericed  evidence  as  to  regularity  of  the  momentum  equation.  Rirther, 
it  is  in  remarkably  close  agreement  with  the  results  obtained  by  Cowley  (1983)  using  a 
series  extension  technique.  In  particular,  Cowley  (1983)  finds  a  singularity  in  the  solution 
at  the  same  time  euad  position  as  the  Lagrangian  computations.  Ingham  (1984)  performed 
an  Eulerian  Fourier  series  expansion  of  the  solution  in  the  direction  along  the  cylinder. 
By  carefully  increasing  the  order  of  expansion  as  the  spectrum  expands  due  to  the  incipi¬ 
ent  singularity,  he  obtained  results  in  close  agreement  with  those  of  both  Van  Dommelen 
&  Shen  (1980a)  and  Cowley  (1983).  The  fact  that  these  three  very  different  procedures 
were  found  to  produce  results  in  excellent  agreement  with  one  another  until  very  close  to 
the  breakdown  of  the  solution  at  separation  is  reassuring,  since  a  number  of  more  con¬ 
ventional  finite  difference  computations  (e.g.  Telionis  &  Tsahalis  (1974),  Wang  (1979), 
Cebeci  (1986))  give  significantly  different  results.  Yet  the  results  of  Henkes  &  Veldman 
(1987)  remain  in  agreement  with  the  three  unronvential  methods  until  relatively  close  to 
the  singularity,  but  disagree  with  Cebeci  (1986)  at  a  signiHcantly  earlier  time.  One  of  the 
difficulties  with  conventional  finite  difference  procedures,  as  pointed  out  by  Cebeci  (1986), 
is  the  need  to  satisfy  the  CFL  condition,  a  condition  which  is  implicitly  satisfied  by  the 
three  procedures  of  Van  Dommelen  &  Shen,  Cowley  and  Ingham. 

Clearly  in  any  numerical  Lagrangian  computations,  it  is  not  possible  to  prove  that  the 
solution  is  regular,  since  the  inevitable  upper  limit  on  resolution  means  that  high  order 
singularities  are  difficult  to  resolve.  However,  in  the  accompanying  numerical  study,  part 
2,  the  boundary  layer  at  the  equatorial  plane  of  a  spinning  sphere  is  solved  using  up  to 
1000  mesh  points  across  the  boundary  layer.  Even  at  such  high  resolution,  no  trace  of 
singular  behavior  wais  observed,  and  derivatives  of  high  order  could  be  evaluated  precisely. 

When  the  fact  that  solutions  to  the  momentum  equations  are  regular  is  accepted, 
(and  for  compressible  flow  in  addition  the  density  must  be  regular),  the  next  question  to 
arise  is  what  implications  such  regularity  has  for  the  structure  of  the  separation  process. 
First,  only  the  continuity  equation  can  develop  singular  behavior,  and  from  (2.2)  or  (2.6) 
it  follows  that  this  is  only  possible  if  the  Lagrangian  gradients  of  x  and  z  become  parallel, 
i.e.  if  at  some  point  s 

Vx  =  A,Vz  ,  (2.7a) 

where  A,  is  a  constant.  Generally,  the  point  s  of  interest  is  the  particle  and  time  at  which 
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(2.7a)  is  satisfied  for  the  first  time.  The  condition  (2.7a)  is  a  three-dimensional  extension 
of  the  two-dimensional  condition  first  pointed  out  by  Van  Dommelen  &  Shen  (1980a);  it 
requires  that  a  Lagrangian  stationary  point,  Vn  =  0,  exists  for  an  oblique  coordinate 

n  =  x-X,z  .  (2.76) 

An  alternate  way  to  phrase  the  condition  for  singular  y  is  to  define  a  unit  vector  in 
n-direction  tangential  to  the  wall, 

in  which  case  a  singularity  occurs  when,  for  all  infinitesimal  changes  in  fluid  particle, 

n.ax  =  0  ,  dx  =  {dx,dz)  .  (2.7d,c) 

This  implies  that  an  infinitesimsj  particle  volume  d^drfd^  around  point  s  has  been  com¬ 
pressed  to  zero  physicsd  size  in  the  n-direction.  But  since  particle  volume  (or  msiss  in 
compressible  flow)  is  conserved,  this  compression  in  the  n-direction  along  the  wall  is  com¬ 
pensated  for  by  a  rapid  expansion  in  the  y-direction  (see  figure  1),  which  drives  the  fluid 
above  the  compressed  region  d^dr}d(  ‘far’  from  the  wall  to  form  a  separating  vorticity 
layer. 

FVom  (2.6)  it  can  be  shown  that  this  process  constitutes  separation  in  the  sense  of  Sears 
&  Telionis  (1975),  since  the  particle  distance  from  the  wall  becomes  too  large,  ‘infinite’,  to 
be  described  on  the  usual  boundary-layer  scale.  Note  that  the  assumed  regulwity  of  x  and 
z  does  not  allow  an  infinite  expansion  in  the  direction  parallel  to  the  wall  but  normsJ  to  n; 
the  particle  can  only  expand  strongly  in  the  direction  away  from  the  wall.  Similarly  for  a 
compressible  fluid,  the  assumed  regularity  of  p  is  inconsistent  with  an  infinite  compression 
of  the  particle  volume.  (At  present  there  is  no  direct  numerical  evidence  for  the  regularity 
assumption  in  the  compressible  case,  although  it  is  of  course  self-consistent). 

From  (2.7)  we  can  derive  generalized  so-called  Moore- Rott-Sears  (MRS)  conditions  at 
the  stationary  point,  similar  to  the  conditions  formulated  by  Sears  &  Telionis  (1975)  for 
two-dimensional  flow.  The  form  of  the  Eulerian  Dy  operator  (2.4)  implies  using  (2.26)  and 
(2.7a)  that  the  vorticity  vanishes  at  that  point,  i.e. 

DyU  =  DyW  =  0  at  Vn  =  0  .  (2.8a) 
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In  fact,  the  Dy -operator  vanishes  for  all  quantities  which  lemsun  regular  in  the  Lagrangian 
domain. 

The  second  Moore- Rott-Sears  condition  is  more  complicated.  Since  (2.7a)  is  equivalent 
to  two  conditions  on  the  Lagrangian  derivatives  of  x  and  z,  in  three  dimensioned  space  we 
expect  it  to  be  satisfied  on  a  curve  of  particles  for  times  beyond  the  first  occurrence  of 
separation  (c.f.  section  3  euid  sub-section  4c).  The  Eulerian  projection  of  the  singular 
curve  on  the  wall  will  be  denoted  by  xmrs  =  {^mrSi^mrs)  a^id  the  Moore-Rott-Sears 
condition  concerns  the  motion  of  this  projected  curve.  To  derive  it,  we  focus  attention  on 
an  arbitrary  point  s  on  the  singular  curve  (rather  than  our  usual  choice  in  which  s  is  the 
first  point  at  which  a  singularity  occurs).  First  we  consider  a  Lagrangian  differential 
along  the  singular  curve  passing  through  point  s,  keeping  time  constant.  Since  x  and  z 
are  functions  of  {  and  t  oidy,  corresponds  to  a  change  in  Eulerian  position  along  the 
projected  curve  which  satisfies  (2.7d), 


n  •  dxMRs  =  0  , 


(2.86) 


so  that  the  singular  curve  is  normed  to  the  local  vector  n.  As  for  any  curve,  the  propagation 
velocity  of  this  curve  is  given  by  the  component  of  the  propagation  velocity  of  points  on 
the  curve  in  the  direction  normal  to  the  curve.  To  find  an  expression  for  it,  we  now 
consider  a  total  differential  in  Lagrangian  space-time  at  the  point  j,  resulting  in  changes 
d®MHS  =  aiid  dzMRS  =  ^^MRS  +  z«d<.  Since  {dx^RSi^^MRs)  satisfies 

(2.8b), 

n - — —  =  n-UMfl5  ,  UAfH5  =  (®»,2»)  ,  (2.8c, d) 

which  shows  that  the  propagation  velocity  of  the  singular  curve  equals  the  flow  velocity  of 
the  singular  particle  a  at  the  considered  position  (®mrs>^mrs)- 

While  this  three-dimensional  form  of  the  MRS  condition  seems  new,  the  general  ap¬ 
plicability  of  the  two-dimensional  case  is  fairly  well  established  both  theoretically  (Moore 
1958,  Sears  &;  Telionis  1975,  Williams  1977,  Shen  1978,  Sychev  1979,  1980,  Van  Domme- 
len  k  Shen  1980b,  1982,  1983a, b.  Van  Dommelen  1981)  and  experimentally  (Ludwig  1964, 
Didden  &  Ho  1985). 

We  can  also  verify  the  notion  of  Sears  &  Telionis  (1975)  that  unsteady  separation 
occurs  in  the  middle  of  the  boundary  layer  rather  than  at  the  wall.  In  the  absence  of  a 
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transpiration  velocity,  the  motion  of  points  on  the  wall  equals  the  motion  of  the  boundary- 
layer  particles  at  the  wall,  cf.  (2.5)  and  (2.3a,  b).  Thus  a  fluid  particle  at  the  wall  can  only 
contract  to  vanishing  size  in  the  n-direction  if  the  wall  itself  performs  the  same  contraction, 
which  is  not  possible  for  a  solid  wall. 

In  the  next  sections  the  nature  of  the  separation  process  is  analyzed.  First  we  form 
local  Taylor  series  expansions  for  the  regular  solutions  to  the  momentum  equations  near 
the  stationary  point,  and  then  we  expand  the  solutions  of  the  continuity  equation  in  an 
asymptotic  series.  This  procedure  is  similar  to  the  one  followed  by  Van  Dommelen  &: 
Shen  (1982)  for  two-dimensional  separation.  In  contrast  to  the  steady  viscous  singularities 
of  Goldstein  (1948)  and  Brown  (1965),  and  the  ideas  of  Sears  &  Telionis  (1975),  the 
unsteady  singularity  is  essentially  inviscid  in  character  and  consists  of  two  vortex  sheets 
separated  by  an  increasingly  large  central  inviscid  region  (as  found  by  Ockendon  (1972) 
for  a  rotating  disc  with  suction,  and  by  Sychev  (1979,  1980),  Van  Dommelen  &  Shen 
(1980b, 1983a,b),  Williams  &  Stewartson  (1983)  and  Elliott,  Cowley  &  Smith  (1983)  for 
steady  separation  over  up-  and  down-stream  moving  walls).  The  leading  order  asymptotic 
structure  of  the  unsteady  singularity  has  also  been  recovered  by  Van  Dommelen  (1981)  as 
a  matched  asymptotic  solution  to  the  Eulerian  boundary-layer  equations.  More  generally, 
Elliott  et  al.  (1983)  showed  that  there  is  a  certain  amount  of  arbitrariness  in  the  Eulerian 
expansions.  The  Lagrangian  expansion  resolves  such  arbitrariness  by  the  assumption, 
(supported  by  various  numerical  data,  see  Van  Dommelen  &  Shen  (1982),  the  closing 
remarks  of  subsection  4c,  and  part  2),  that  the  leading  order  coefficients  in  the  Taylor 
series  expansion  near  the  stationary  point  are  non-zero. 

3.  Three-dimensional  separation  singularities 

In  this  section  we  find  the  leading  order  term  of  an  asymptotic  analysis  which  describes 
the  local  structure  of  the  flow  cS  unsteady  separation  is  approached.  The  time  and  position 
at  which  the  separation  singularity  first  develops  will  be  denoted  by  the  subscript  a,  thus 
for  example 

(Vn),  =  0  ,  (3.1a) 

where  n  is  the  oblique  coordinate  corresponding  to  the  initial  separation,  defined  in  (2.76) 
as 
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n  =  X  —  \gZ 


(3.U) 


Note  that  the  definition  of  the  x-  and  z-coordinates  can  simply  be  interchanged  if  n  and 
z  are  not  independent  coordinates.  In  index  notation,  (3.1o)  can  be  written  aa  =  0, 
where  we  will  adopt  the  convention  to  omit  the  subscripts  comma  (to  indicate  Lagrangian 
derivatives)  and  a  (to  indicate  the  separation  particle  at  the  separation  time)  if  they  occur 
together  (i.e.  rn  =  (ni)J. 

The  solution  of  the  continuity  equation  (2.2)  for  y  can  be  greatly  simplified  by  a 
number  of  coordinate  transformations  for  both  the  particle  position  coordinates  (®,«)  and 
the  Lagrangian  coordinates  (^,^7,  C)-  Here  we  will  select  transformations  which  preserve  the 
Jacobian  J  (2.26),  since  these  are  algebraically  more  simple  than  transformations  which 
preserve  the  physical  volume  HJ,  or  mass  pHJ. 

As  a  first  trsmsformation,  we  drop  the  position  coordinate  x  in  favor  of  n,  shift  the 
Lagrangian  coordinate  system  to  the  separation  particle  s,  and  rotate  it,  resulting  in  the 
set  of  coordinates 


3 

n  —  X  —  \gZ  ,  z  ,  hi  —  <iij (^j  ~  ^j$)  )  (3.2(1, 6, c) 


where  Oij  is  an  orthonormal  rotation  matrix  which  is  is  chosen  to  eliminate  the  mixed 
derivatives  nu,  nja,  and  7x23.  Therefore,  expanding  n  and  2  in  a  Taylor  series  expansion 
about  the  separation  point,  we  obtain 

+  .. .  +  + 

i=l  ' 


+ 


(3.3a) 


3 

2  =  2, -t- ^  2i/c, •  + .. .  +  6<2,  4- . . .  ,  (3.36) 

t=i 

where  8t  ■=  t  —  tg. 

However,  if  tg  is  the  first  time  that  a  stationary  point  occurs,  the  Taylor  series  co¬ 
efficients  in  (3.3)  cannot  be  completely  arbitrary:  the  singularity  condition  may  not  be 
satisfied  anywhere  for  6t  <  0.  The  condition  for  a  singularity  to  exist  for  earlier  times  at 
some  neighboring  point  is,  in  terms  of  n  and  2, 


)^>»  —  ® 


(3.4a) 
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where  A  is  the  ratio  between  Vz  and  V®  at  the  neighboring  singular  point.  Elxpanding 
(3.4a)  in  a  Thylor  series,  we  obtain 


fHiki  —  ZiSX  +  riiSt  +  . . .  =  Q  for  1  =  1,2, 3  ,  (3.46) 

where  6A  =  A  —  A,.  If  all  three  coefficients  nn,  n2j,  and  7133  were  non-zero,  a  solution 
to  (3.46)  could  be  found  for  St  <  0,  contradicting  our  assumption  that  the  singularity 
develops  first  at  St  =  0.  (Strictly,  because  of  the  higher  order  terms  omitted  in  (3.46), 
the  solution  must  be  found  iteratively,  however,  the  iterations  converge  because  the  higher 
order  terms  act  as  a  contraction  mapping  sufficiently  close  to  point  s).  Therefore  at  the 
first  occurrence  of  separation,  at  least  one  of  7111,7132,  or  7133  must  be  zero,  and  we  will 
reorder  (Ari , A:2 , Ajs )  such  that  Tin  vanishes.  In  euidition,  the  coefficients  7133,7133,21  cannot 
all  be  non-zero,  since  by  solving  for  6A,  and  ki,  it  again  follows  that  a  angtilarity  exists 
for  St  <  0.  Without  loss  of  generality,  we  assume  that  zi  is  zero,  since  if  either  7133  or  7133 
vanishes,  the  (ki^kj^ks)  coordinate  system  can  be  rotated  further  to  eliminate  zi. 

It  follows  that  in  some  suitably  oriented  Lagrangian  coordinate  system  the  conditions 
Till  =  =0  8^"®  necessary  at  the  time  when  separation  starts.  This  implies  two  addi- 

tionaJ  conditions  on  x(^,t)  and  z{^,t),  besides  the  two  conditions  implicit  in  (3.1a).  Since 
Lagrangian  space-time  is  four-dimensional,  in  general  we  do  not  expect  that  more  than 
four  conditions  can  be  satisfied  at  any  time.  Hence,  in  the  remainder  of  this  section  we 
will  assume  that  the  values  of  the  remsuning  derivatives  can  be  completely  arbitrary  and 
in  general  non-zero. 

However,  when  the  functions  x  and  z  are  not  arbitrary,  but  restricted  by  constraints  of 
symmetry  in  the  flow,  the  latter  assumption  needs  to  be  reconsidered,  since  the  symmetry 
requires  that  various  derivatives  must  vanish.  Examples  are  two-dimensional  flow,  and  the 
flows  discussed  in  the  next  section. 

Under  the  assumption  that  the  remaining  coefficients  in  the  Taylor  series  have  arbi¬ 
trary  values,  the  transformation 


z  =  2-z(^,,0  ,  h  =  n  -  n{i,,t)  - 


ki  =  ki 


ki  = 


^22^3^2  —  ^332:2^3 
y/nl^zj  +nljz^ 


ki  = 


-1-712223^:3 
\/«22^3  +«L^2 


(3.5a, 6) 
(3.5c,  d,e) 


12 


where 


^(2) _ ^3^33 

'  2(712223  +^13322)  ’ 

eliminates  the  7133  term.  The  final  coordinate  transform 


(3.5/) 


/i  =  *1  +  ,  ii  =  fa  ,  <3  =  fa  , 

Tim 


n  =  n  —  —  H,Stz  ,  z  =  z  , 


(3.6a,  6,  c) 


(3.6d,  e) 


where 


^f3)  fijjjn333  —  3niiinii3ni33  +  2njj3 

'''  “  6n=*  53  ’ 

D71jjj23 

eliminates  the  nus,  71333,  and  713  derivatives. 


H,  = 


Tlm7l3  —  71113711 

fill!  23 


(3.6/, j) 


The  transformed  position  coordinate  ti  corresponds  to  an  oblique  coordinate  measured 
from  a  moving,  curved  line  through  the  separation  particle,  viz. 


71  =  X-Xo(2,<)  , 


(3.7a) 


where 


®  =  x-x(^,,<)  ,  z  =  z-z{i„t)  , 

xo(2,<)  =  A,z  +  +  Hs8tz  . 


(3.76,  c) 
(3.7d) 


Note  that  the  curved  line  x  =  xo(2,<),  which  can  be  viewed  as  the  line  along  which 
the  separation  initially  develops  (see  below),  does  not  have  a  singular  shape  at  the  first 
occurrence  of  separation. 

The  Taylor  series  expansions  for  n  and  2  near  the  separation  point  become 
n  =  ^n22f2  +  ^  InijkUljlk  +  . . .  +  ^  nJi  +  .. .  ,  (nna  =  n33-,  =  713  =  0),  (3.8a) 


2  —  22/2  +  23/3  4-  .  •  • 


(3.86) 


The  characteristics  of  the  Jacobian  equation  (2.2)  for  y  are,  in  terms  of  the  new  coordinates, 


^  _  pH 

dy  po  Hq 


{-Z3n22h  +  .  •  •  }  ) 


(3.9a) 
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(3.96) 


^  +  \n133ll  +ni6i)  +  ...}  , 

-^  =  -^-jj-{—Z2{\niiil\  +  ^riiisll  +  riiSt)  +  . . .}  ,  (3.9c) 

with  a  singularity  occurring  when  all  three  right  hand  side  expressions  vanish  (note  that 
not  all  three  are  independent).  Near  the  point  s,  (3.9o)  is  zero  on  a  surface  approximating 
the  I2  =  0  plane,  while  both  (3.96)  and  (3.9c)  vanish  at  points  depending  on  the  nature 
of  the  quadratic  expression  (jnm/J  +  iniaa/^).  If  this  quadratic  is  hyperbolic,  singular 
particles  occur  along  hyperbolic  lines  regardless  of  the  agn  of  6t.  Thus,  if  =  0  is  to  be 
the  first  time  that  separation  occurs,  the  quadratic  must  be  elliptic,  and  of  the  same  sign 
as  the  constant  term  when  6t  <  0.  This  requires  >  0  and  nmni  <  0;  we  will 

choose  the  positive  Zj -direction  such  that 

«22nni  >0  ,  njjfiiaa  >0  ,  njjfii  <0  .  (3.10o,  6,c) 


The  Lagrangian  description  of  the  separation  process  can  now  be  completed  by  the 
determination  of  y  at  times  shortly  before  the  initiaJ  occurrence  of  separation.  At  t  =  t, 
the  boundary-layer  approximation  is  obviously  no  longer  valid  because  from  the  integral 
(2.6)  it  follows  that  y  becomes  infinite  at  the  stationary  point.  However,  the  rate  of  growth 
near  this  point  can  be  found  by  means  of  an  asymptotic  expansion.  To  find  local  scalings, 
we  follow  the  guiding  principles  of  Vein  Dyke  (1975).  In  general,  we  attempt  to  scale  the 
Lagrangian  coordinates  li  and  the  position  coordinates  n,  z  and  y  to  variables  Li,N,Z, 
and  Y  such  that  in  the  inner  region  the  Jacobian  equation  for  Y,  i.e. 

Yl,  n,  Yl, 

Zci  Zlj  Zl, 

has  non-singular  leading  order  coefficients.  This  suggests  that  the  6t  term  in  (3.96),  which 
ensures  the  absence  of  singular  points  for  6t  <  0,  should  be  retained.  Further,  for  =  0 
we  want  to  match  the  solution  close  to  the  stationary  particle  to  a  solution  for  y  which 
is  regular  away  from  this  point.  Thus  we  want  to  retain  those  terms  which  ensure  the 
absence  of  singular  points  away  from  particle  at  time  St  =  0,  i.e.  the  Zj-  and  Zj -terms 
in  (3.96)  and  the  Za-term  in  (3.9a).  The  appropriate  scaling  is  therefore 

Zi=|6Z|U,  ,  l2  =  \St\H2  ,  Z3  =  |6Z|>T3  ,  (3.12a,6,c) 


Po 

pH 


(3.11) 


JciN,Y,Z)  = 
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n  =  |^<|^iV  ,  z  =  \6i\^^Z  ,  y=\6t\-^Y  .  (3.12d,e,/) 

These  scalings  suggest  that  the  separation  process  occurs  in  a  relatively  thin  strip,  n  ~ 
along  a  segment  of  the  separation  line  x  =  ®o(^>0  of  length  z  ~  |^<|  >. 

For  the  scaling  (3.12),  the  solution  for  Y  is  most  easily  found  by  integration  of  (3.11) 
as  in  (3.9a),  where  ij  and  L3  are  eliminated  in  favor  of  N  and  Z,  which  are  constant  along 
the  lines  of  integration,  using  (3.8).  The  result  is: 

P.«.  V-oo  v'P{£i  N,Z)  Jl,  ^P{L-,  N,  Z) 

where 

P(I;  N,Z)  =  -  ifijj  {zlnniL^  +  (3n,33  Z’  -  6^1  z^)I  -  JV}  ,  (3.136) 

and  Lo{N,  Z)  is  the  real  root  of  the  cubic  P  .  This  root  is  a  unique  and  continuous  function 
of  N  and  Z  since  P  is  a  monotonically  decreasing  function  of  L  from  (3.10). 

The  choice  of  agn  of  the  square-root  in  (3.13a),  and  the  limits  of  integration  are 
determined  by  the  topology  of  the  lines  of  constant  N  and  Z.  In  physical  space  these  lines 
are  straight  and  pass  vertically  through  the  boundary  layer;  however,  in  Lagrangian  space 
they  are  highly  curved  near  the  separation  particle,  as  shown  qualitatively  in  Rgure  2a.  The 
lines  can  be  divided  into  three  segments  corresponding  to  three  asymptotic  regions.  The 
lower  segments  start  at  the  wall  and  extend  upward  towards  the  vicinity  of  the  separation 
particle.  Because  the  Jacobian  is  nowhere  singular  along  these  segments,  the  y-positions 
of  the  fluid  particles  remain  finite  on  the  boundary-layer  scale,  i.e.  the  scaled  coordinate 
Y  is  small.  Hence,  these  lower  segments  give  rise  to  a  layer  of  particles  at  the  wall  with 
a  thickness  comparable  to  that  of  the  original  boundary  layer,  this  is  shown  schematically 
in  figure  1. 

Along  the  central  segments,  the  lines  of  constant  N  and  Z  pass  through  the  vicinity 
of  the  separation  particle.  Here  the  y-position  of  the  particles  grows  rapidly,  and  is  given 
in  scaled  form  by  (3.13).  Thus  the  central  segments  give  rise  to  the  intermediate,  thicker, 
layer  of  particles  shown  in  figure  1.  The  topology  of  the  central  segments  in  the  Lagrangian 
domain,  figure  2a,  determines  the  choice  of  sign  in  (3.13a).  From  (3.8)  and  (3.9)  it  follows 
that  on  integrating  upwards,  Li  increases  from  large  negative  values  towards  Lo{N^Z). 


(3.13a) 
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Since  Y  is  increasing  along  this  part,  the  negative  sign  in  (3.13a)  applies.  At  position  Lq, 
the  lines  of  constant  N  and  Z  turn  around  in  the  Lagrangian  domain  and  Lj  again  tends 
to  — oo;  along  this  second  part  the  positive  sign  in  (3.13a)  applies. 

Along  the  third  segments,  the  lines  of  constant  N  and  Z  proceed  upwards  toward  the 
external  flow.  As  in  the  lower  segments,  the  Jacobian  is  no  longer  small  here.  Thus  the 
changes  in  y  are  finite  on  boundary- layer  scale,  and  the  third  segments  give  rise  to  a  layer 
of  particles  with  a  boundary-layer  scale  thickness,  atop  the  central  region,  as  shown  in 
figure  1. 

Hence,  the  separation  structure  is  one  in  which  the  boundary  layer  divides  into  a 
central  layer  of  physical  thickness  proportional  to  between  two  ‘sandwich’ 

layers  of  thickness  proportional  to 

The  structure  (3.13)  is  identical  to  the  one  obtained  by  Van  Dommelen  &  Shen  (1982) 
for  two-dimensional  separation,  except  that  the  coefficients  now  depend  on  the  position  Z 
along  the  describing  line  ico.  A  convenient  way  to  illustrate  the  influence  of  the  position 
Z  is  to  scale  out  the  coefficients  using  a  procedure  similar  to  Van  Dommelen  (1981): 


=/3jii  =^2(Z^  +  l)ilI  ,  L2=  0002^12  =0002^  {z^ ,  (3.14a,6) 

N  =  a0lN  =  a0l{Z'^  -f  l)>Ar'  ,  T  =  ^  =  — - -  ,  Z  =  ^Z  ,  (3.14c,d,e) 

where  the  tilde- variables  scale  out  the  Taylor  series  coefficients,  the  starred  variables  scale 
out  Z,  and 


a  =  ,  1  =  {\z\'^22nu\) 


1  p,Hs 

PO,ffos 


00  = 


23^111 


,  /31  = 


In  terms  of  the  starred  variables  (3.13)  reduces  to 


^111  / 


^111 


(3.14d,e) 
(3.14/,  5,  h) 


r^o  dL*  ^  _ dL* 

1-00  -  31*  -  T*’  ^27V*  -  3£*  -  L* 


(3.14t) 


L;(iV*)  =  f(iV*)  , 


where 
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(3.15a) 


and  the  function  I  is  the  inverse  to  the  cubic  N*  =  5/®  +  i.e. 

J(iV)  =  (jv* +(1 +  +  (jv* -(1  +  .  (3.156) 

The  values  of  and  depend  on  the  choice  of  the  Eulerian  coordinates  (*,«), 

but  not  on  the  definition  of  the  Lagrangian  coordinates. 


An  alternative  expression  for  Y*  can  be  found  in  terms  of  the  incomplete  elliptic 
integral  of  the  first  kind  F{<f>\m): 


(3.I60) 


where 


A(Ar-)  =  (3(£;’  +  1))  * 


(p{Ll^N*)  =  2  arctan 


tMN  )  =  -  H - ?■ 

^  ’  2  4A2 


(3.166,  c) 
(3.16ci) 


Elliptic  integrtds  are  distorted  identity  functions,  (in  particulsir  F(^|0)  =  ip  exactly),  so 
that  the  arctan  is  responsible  for  the  major  variations  in  Y  along  the  characteristics. 


Farther  terms  in  the  asymptotic  expansions  (3.8)  and  (3.16)  can  be  found  in  princi¬ 
ple.  We  note  that  the  next  term  in  the  expression  for  Y  does  not  involve  a  logarithmic 
correction,  even  though  logarithmic  second  order  terms  do  arise  for  the  symmetric  flows 
studied  in  the  next  section. 


We  now  turn  to  the  physical  interpretation  of  these  results.  The  boundary-layer 
thickness  is  asymptotically  determined  by  the  position  of  the  upper  particle  layer  in  figure 
1;  letting  LJ  —*  —00  along  the  positive  branch  of  (3.16a),  we  obtain  the  scaled  boundary- 
layer  thickness  as 


(3.17) 


The  function  Y'^*{N*)  gives  the  general  shape  of  the  boundary-layer  thickness  in  a  cross- 
section  of  constant  z.  For  large  values  of  N*  the  boundary-layer  thickness  decays  toward 
zero  much  more  slowly  than  suggested  by  the  sketch  in  figure  1.  Nevertheless,  at  the  outer 
edges  of  the  thin  separation  region,  the  solution  still  matches  with  finite  values  of  y;  for 
from  (3.12),  (3.14),  and  (3.16) 

,  4a  i 


y 


4aS  1  ,  y/3\  1  ,  , 


(3.18) 
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To  show  the  dependence  of  the  boundary-layer  thickness  on  the  coordinate  z,  contours 
of  constant  y"*"  in  the  7V,Z-plane  axe  plotted  in  figure  2b.  Note  that  the  coordinate  N 
is  measured  from  the  oblique,  curved,  separation  line.  Actual  lines  of  constant  boundary- 
layer  thickness  might,  for  example,  appear  as  sketched  in  figure  2c,  which  has  been  drawn 
by  taking  |5<|  =  0.06  and  unit  values  for  various  coefficients  in  (3.7)  and  (3.14).  Asymptot¬ 
ically,  the  boundary-layer  thickness  has  the  shape  of  a  crescent  shaped  ridge.  The  crescent 
shape  is  long  and  thin,  i.e.  quasi-two-dimensional,  because  from  (3.12)  the  n  length  scale 
is  asymptotically  shorter  than  the  z  length  scale  (note  that  for  three-dimensional  steady 
separation  Smith  (1978)  has  proposed  a  quasi-two-dimensional  structure).  In  an  Eulerian 
numerical  calculation,  the  development  of  such  a  crescent-shaped  ridge  may  be  a  possible 
diagnostic  indicating  the  presence  of  a  singularity. 

Evidence  of  this  type  of  singularity  is  provided  by  Ragab’s  (1986)  calculations  for 
impulsively  started  flow  past  a  4:1  prolate  spheroid  inclined  at  a  30“  angle  of  attack.  His 
results  strongly  suggest  that  the  displacement  thickness  becomes  unbounded  away  from 
the  symmetry  line.  However,  it  is  not  possible  to  deduce  the  shape  of  the  singularity  from 
the  results  presented. 

A  point  of  interest  is  the  decay  of  the  boundary-layer  thickness  along  the  describing 
line  for  large  Z.  FVom  (3.12),  (3.14),  and  (3.16), 


for  |^<|  >  <  |z|  <  1 


(3.19) 


Hence  for  increasing  z,  the  separation  structure  expands  in  n-direction,  while  the  thickness 
of  the  boundary-layer  decreases. 

The  particle  propagation  velocity  n  which  ^ves  rise  to  the  accumulation  of  particles 
at  the  separation  line  is,  according  to  (3.8a),  given  to  leading  order  by 


n 


|(J<|  ’n,Tj 


(3.20) 


To  describe  this  in  the  more  familiar  Eulerian  coordinates,  the  transcendental  relationship 
(3.16)  must  be  inverted  to  the  form 


Ll=L]{N%Y-) 


(3.21) 
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The  inversion  has  been  performed  numerically,  and  in  figure  2d  we  present  contours  of  LJ 
in  the  plane.  fVom  (3.14), 


+  .  (3.22) 

It  follows  that  the  lines  of  constant  shown  in  figure  2d  describe  the  shape  of  the  lines 
of  constant  E  in  cross-sections  of  constant  z  through  the  separation  structure.  They  also 
give  the  asymptotic  shape  of  the  lines  of  constant  velocity  components  x  and  i  euid  density 
p  in  these  cross-sections,  since 

{x,z,p)  =  (x,,i,,p,)  +  \6t\^  (^{xuZuPi)^2{Z^  +  l)^L^  +  {x3,Z3,P3)-^-Z^  4-  ...  . 

(3.23) 

We  note  that  the  topology  of  figure  2d  for  |^<|  0  seems  quite  close  to  the  computed  lines 

of  constant  velocity  presented  by  Van  Domm<;len  (1981)  for  finite  |5<|,  and  thus  lines  of 
constant  velocity  might  be  a  useful  indication  of  an  incipient  unsteady  separation. 

The  next  point  of  interest  is  the  shape  of  the  velocity  profiles.  According  to  (3.23), 
in  Eulerian  space  the  velocity  profiles  must  develop  a  large  flat  region  of  nearly  constant 
velocity  as  separation  is  approached.  However,  accepting  the  numerical  results  of  Van 
Dommelen  (1981),  this  flat  region  is  only  evident  extremely  close  to  the  singularity,  so  that 
resolution  problems  or  finite  Re’  lolds  number  effects  tend  to  obscure  the  phenomenon. 
From  (3.23)  and  figure  2d,  the  velocity  profiles  near  an  incipient  three-dimensional  separa¬ 
tion  must  have  a  local  maximum  or  minimum  in  velocity.  However,  this  is  not  necessarily  a 
precise  indication  of  incipient  separation.  For  example,  in  the  case  of  the  circular  cylinder, 
a  minimum  in  the  velocity  profiles  develops  relatively  quickly,  after  |  diameter  motion, 
yet  separation  occurs  much  later,  after  |  diameter  motion.  Figure  2e  shows  the  shape  of 
the  velocity  profiles  near  the  interior  extrema.  The  shapes  of  the  velocity  profiles  in  the 
sandwich  layers  at  the  edges  of  figure  2e  cannot  be  found  from  asymptotic  analysis  since 
they  depend  on  the  precise  details  of  the  earlier  evolution  (cf.  the  remarks  below  (3.24) 
and  part  2). 

A  more  significant  sign  of  the  start  of  separation  might  be  a  transverse  expansion 
of  the  lines  of  constant  vorticity  near  the  velocity  minimum/maximum;  since  the  above 
analysis  is  inviscid  to  leading  order,  the  vorticity  lines  closely  follow  the  motion  of  the 
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boundary-layer  particles.  In  the  boundary-layer  approximation,  the  vorticity  is  the  y- 
derivative  of  the  velocity  distribution.  The  corresponding  asymptotic  topology  of  contours 
of  dL\ldY*  is  shown  in  figure  2f.  This  topology  seems  close  to  the  computed  vorticity 
Hues  presented  by  Van  Dommelen  (1981)  for  a  time  near  separation. 

The  asymptotic  structures  of  the  upper  and  lower  vorticity  layers  are  similar  to  the 
two-dimensional  case  (Van  Dommelen  1981).  Expressed  in  terms  of  Eulerian  coordinates, 
they  take  the  form  of  regular  Taylor  expansions: 


mnr>0 


and 


(3.24a) 


Y.  .  (3.246) 

mnr>0 

respectively,  where  the  sums  run  over  the  non-negative  integers,  and  the  Prandtl  transfor¬ 
mation  y  —  y  —  y'^ {x ,2,6t),  describes  the  motion  of  the  upper  layer. 

Substituting  (3.24)  into  the  boundary-layer  equations,  we  find  that  the 
Pmnr>  (»Ti,n  >  0,  r  >  1)  and  the  (Tn,n,r  >  0)  are  determined  in  terms  of  the 

(^mno>^mnotPmno)>  these  latter  functions  are  indeterminate  due  to  the  depen¬ 

dence  of  the  solution  on  earlier  times.  The  (^^mno’^mno’^mno)  "^^st,  however,  satisfy  the 
boundary  conditions  (2.5a)  at  the  wall,  and  match  both  at  the  outer  edge  of  the  boundary 
layer  (see  (2.56)),  and  with  the  central  inviscid  low-vorticity  region.  At  fixed  N  and  Z, 
the  latter  matching  conditions  yield  from  inverting  (3.16)  and  using  (3.23), 


(^000  ’  ^000  ^Pooo)  Pa)  Pi)  y  +00  ,  (3.25a) 

7^  r 

(“^oo>^o^oo’/0o^oo)  (=^3,2,,^,)  -  (^^i,2i,Pi ^  as  y  -  -oo  .  (3.256) 

7  y 

Asymptotic  matching  conditions  can  also  be  derived  os  |A'|,|Z1  — *  oo,  as  Van  Dommelen 
(1981)  has  done  for  two-dimensional  flows. 


A  final  point  of  interest  is  the  ‘accessibility'  of  the  region  of  flow  beyond  the  time 
of  initial  separation.  In  a  steady  Eulerian  computation,  Cebeci,  Khattab  &  Stewartson 
(1981)  took  the  accessible  region  to  be  the  domain  where  a  boundary-layer  solution  can  be 
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found,  (whether  it  is  still  an  asymptotically  correct  solution  of  the  Navier-Stokes  equations 
in  the  presence  of  interaction  or  not).  In  the  Lagrangian  case,  some  care  is  needed,  because 
the  singular  continuity  equation  is  integrated  separately.  Numerical  experiments  such  as 
the  one  in  part  2  do  in  fact  suggest  that  the  non-singular  momentum  equations  can  be 
integrated  pwt  the  separation  singularity  without  apparent  difficulty.  When  that  is  done, 
the  vertical  lines  through  the  boundary  layer  appear  in  Lagrangian  space  as  shown  in  figure 
2g  rather  than  figure  2a.  For  the  shaded  particles,  y  is  indeterminate;  these  particles  may 
be  thought  of  as  having  disappeared  at  infinite  y.  Yet  the  continuity  equation  can  still 
be  integrated  along  all  lines  of  constant  n  and  z  which  start  at  the  wall.  A  singularity 
develops  only  on  the  line  passing  through  the  saddle  point  in  figure  2g,  which  for  0  <  •C  1 
corresponds  to  a  singular  line  segment 

iV  -  -f-  (l  -  j  *  .  (3.26) 

However,  the  solution  so  obtained  must  be  considered  meaningless  at  least  for  all  particles 
which  have  at  some  previous  time  passed  through  the:  singular  curve.  For  that  reason,  we 
define  the  region  of  inaccessibility  as  those  stations  {x,z)  which  contain  particles  which 
have  at  any  time  been  on  the  singular  curve.  Initially,  the  region  of  inaccessibility  will 
primarily  expand  in  the  z-direction  through  the  scaling  (3.12e).  In  the  n-direction  it  will 
expand  by  means  of  the  motion  of  the  describing  line  (3.7)  and  additionally  through  the 
motion  of  the  particles  which  propagate  downstream  away  from  the  singular  curve.  Thus, 
the  region  of  inaccessibility  extends  over  a  finite  surface  area,  rather  than  just  the  curve 
(3.26),  in  agreement  with  the  steady  Eulerian  definition  of  Cebeci  et  al.  (1981). 

Naturally,  the  singularity  structure  derived  here  will  not  remain  asymptotically  correct 
arbitrarily  close  to  <  =  <, ,  because  the  normal  velocity  above  the  central  inviscid  region 
becomes  infinite  at  t  =  t, .  From  a  study  of  the  Navier-Stokes  equations  it  is  found  that  the 
singularity  is  smoothed  out  when  a  ‘triple-deck’  interaction  comes  into  operation  for  bt  — 
0{Re~^),  at  which  point  the  scaled  boundary-layer  thickness  is  C>(/?e”).  Because  the 
singularity  is  quasi-two-dimensional,  the  scalings  and  governing  equations  are  essentially 
those  derived  by  Elliott  et  al.  (1983)  for  two-dimensional  flows,  but  with  the  addition  of  a 
passive  z-momentum  equation.  In  the  central  interaction  problem,  the  coordinate  z,  which 
has  an  interaction  length  scale  0{Re~  ^  ),  only  appears  as  a  parameter.  However,  it  is  not 
dear  whether  the  singularity  will  be  completely  removed  by  the  interaction  (Smith  1987). 
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4.  Three-dimensional  symmetric  separation 

In  the  previous  section  a  separation  singularity  structure  wm  derived  assuming  that 
the  flow  was  arbitrary,  an  assumption  that  might  be  appropriate  for  flow  past  an  asymmet¬ 
ric  body.  However,  in  the  case  of  a  spheroid  at  relatively  small  angles  of  attack  it  is  likely 
that  separation  first  occurs  on  one  of  the  symmetry  lines;  indeed  numerical  calculations 
confined  to  the  symmetry  line  have  been  performed  on  this  basis  (Wang  &  Fan  (1982) 
Cebeci,  Stewwtson  &  Schimke  (1984)).  In  sub-section  (a)  below  we  derive  the  form  of  the 
singulzuity  appropriate  for  separating  flows  where  the  separation  line  crosses  a  symmetry 
line  normally. 

However,  this  is  not  the  only  type  of  symmetric  separation  of  interest.  When  a  sphere 
is  impulsively  rotated  about  a  diameter,  centripetal  effects  generate  a  boundajy-layer  flow 
towards  the  equator.  After  a  finite  time  an  equatorial  singularity  develops  as  a  result  of 
a  boundary-layer  collision.  The  structure  of  this  angularity  on  the  symmetry  line  has 
been  determined  by  Banks  &  Zaturska  (1979),  and  Simpson  &  Stewartson  (1982a).  In 
this  case  the  sepzu’ation  line  coincides  with  the  symmetry  line.  Similar  singulairities  occur 
after  a  finite  time  at  the  apex  of  a  horizontal  circular  cylinder  which  is  impulsively  heated 
(Simpson  &:  Stewartson  1982b),  at  the  inner  bend  of  a  uniformly  curved  pipe  through  which 
flow  is  impulsively  started  (Lam  1988),  and  at  the  stagnation  points  on  a  two-dimensional 
cylinder  in  oscillating  flow  as  a  result  of  steady  streaming  effects  (Vasantha  &  Riley  1988). 

A  more  general  form  of  the  singularity  generated  by  two  symmetric  colliding  boundary 
layers  on  a  smooth  wall  would  first  develop  at  a  point  rather  than  along  the  entire  symmetry 
line.  For  example,  such  a  singularity  might  develop  on  the  equator  of  an  ellipsoid  which 
is  rotated  about  one  of  its  principal  axes,  or  in  starting  flow  through  a  curved  pipe  with 
non-uniform  curvature,  or  at  the  apex  of  a  heated  ellipsoid.  In  sub-section  (b)  the  three- 
dimensional  structure  of  such  a  singularity  is  derived.  The  results  on  the  symmetry  line 
agree  with  previous  authors,  but  the  simplicity  of  the  Lagrangian  approach  allows  us  to 
determine  additionally  the  singularity  structure  off  this  line.  The  latter  is  a  necessary 
preliminary  in  order  to  formulate  subsequent  asymptotic  stages  in  the  separation  process. 

Another  class  of  separation  singularities  are  rotationally  symmetric  about  the  sepa¬ 
ration  point,  so  that  the  separation  line  degenerates  to  a  point.  For  example,  singularities 
develop  after  a  finite  time  on  the  axis  of  a  spinning  disc  or  sphere  whose  direction  of 
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rotation  is  impulsively  reversed  (Bodonyi  &  Stewartson  1977,  Banks  &  Zaturska  1981, 
Stewartson,  Simpson  &  Bodonyi  1982,  Van  Dommelen  1987),  and  at  the  apex  of  a  sphere 
which  is  impulsively  heated  (Brown  &  Simpson  1982,  Awang  &  Riley  1983).  The  structures 
of  these  singularities,  which  differ  due  to  the  presence  and  absence  of  swirl,  are  derived  in 
sub-sections  (c)  and  (d)  respectively.  The  results  on  the  axis  agree  with  those  of  previous 
authors,  while  the  singularity  structures  off  the  axis  are  new. 


(a)  Lateral  symmetry 

When  the  boundary-layer  flow  is  symmetrical  about  a  line  along  the  surface  of  the 
body,  the  describing  line  of  separation  must  either  cross  the  symmetry  line  normally  or 
coincide  with  it.  In  this  sub— section  we  will  address  the  case  of  normal  crossing,  leaving 
the  second  possibility  to  the  next  sub-section. 

For  consistency  with  section  3,  we  identify  the  compressed  coordinate  n  with  *  and 
take  the  77-plane  as  the  symmetry  plane  so  that  x  is  an  even  function  of  ^  and  z  an  odd 
function.  Then  the  analysis  is  a  simpler  version  of  the  one  in  the  previous  section.  The 
only  transformation  of  the  Lagrangian  coordinate  system  needed  is  a  rotation  around  the 
(^-axis  to  eliminate  the  xjj  derivative.  Also,  the  discussion  concerning  which  derivatives 
must  be  zero  if  is  the  first  separation  time  (see  (3.4)  and  following)  can  be  restricted  to 
the  symmetry  plane  to  show  that  the  second  order  derivative  which  is  forced  to  be  zero 
must  lie  within  the  symmetry  plane. 


Hence  the  structure  of  the  separation  process  remains  basically  unchanged,  although 
the  describing  line  of  separation  simplifies,  and  is  now  symmetric  about  the  symmetry  line 
z  =  0  (cf.  (3.7)): 

n  =  X  -  x(^,,<)  -  ^2^  .  (4.1) 


A  degenerate  case  is  two-dimensional  flow,  where  x  is  totally  independent  of  (,  and  the 
separation  line  becomes  a  straight  generator  in  the  z-direction.  In  addition,  the  coefficient 
/?!  vanishes,  which  suppresses  the  decay  of  the  boundary-layer  thickness  with  2.  The 
resulting  structure  is  described  in  detail  by  Van  Dommelen  (1981). 

Thus  lateral  symmetry,  or  more  strongly  two-dimensionality,  does  not  fundamentally 
alter  the  separation  process.  This  conclusion  is  consistent  with  the  symmetry  line  calcula¬ 
tions  of  Cebeci,  Stewartson  &  Schimke  (1984). 
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(b)  Symmetric  boundary-layer  collision 

When  the  describing  line  coincides  with  the  symmetry  line,  signiflcant  changes  in 
structure  are  unavoidable,  since  the  flow  is  symmetric  while  the  separation  structure  illus¬ 
trated  in  figure  2  is  asymmetrical. 

We  identify  the  compressed  coordinate  n  again  with  x,  but  now  we  assume  that  the 
7/,  ^-plane  is  the  symmetry  plane,  so  that  a;  is  an  odd  function  of  (  while  z  is  an  even 
function.  A  singularity  occurs  when  first  vanishes  at  the  symmetry  plane,  since  the 
derivatives  x^^  and  are  zero  by  symmetry.  Since  the  first  occurrence  of  a  zero  value 
must  occur  where  x^^  is  a  minimum,  the  second  order  derivatives  x^^  and  x^^  must  vanish, 
while  the  other  second  order  derivatives  are  zero  by  symmetry. 

I'he  fact  that  all  the  second  order  derivatives  are  zero  invalidates  the  scalings  for 
77  and  y  made  in  the  previous  section  (e.g.  (3.12),  (3.14)),  hence  a  separate  analysis 
with  signiflcant  modifications  is  needed.  Proceeding  along  similar  lines  as  in  the  previous 
section,  a  local  Lagrangian  coordinate  system  ki,k2,k3  is  introduced  with  OTigin  at  the 
separation  particle,  but  with  the  same  orientation  as  the  original  axis  system.  A  rotation 
of  this  coordinate  system  around  the  fej-axis. 


ki  =  (  ,  ^2  = 


-*3^2  —  Z2k3 


^2^2  +  ^3^3 


'4  +  4 


x  =  x  ,  z  =  z-z($,,<)  , 


can  be  made  to  eliminate  the  zj -derivative.  The  shearing  transformation 


^  ^  1  I2  —  k2  +  - — ^^3 

®122 


/3  =  k. 


x  =  x  ,  z  =  z-z{^,,t)  , 

eliminates  the  X123  derivative,  resulting  in  the  Taylor  series  expansions 


rXii  j  /]  4-  +  2^i33^\^l  +  •  •  •  +  ^tx  ]  /j  + 


^3^3  +  •  •  • 


(4.2a,  6,  c) 


(4.2d,e) 


(4.3a, b, c) 

(4.3d,  e) 


(4.4a) 


The  expressions  for  the  characteristics  of  the  Jacobian  equation  for  y  become 


^  _  pff 

dy  po  Ho 


U  2:3X122^2  -!-•  ••}  , 


(4.46) 


(4.5a) 
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(4.56) 


^  ~  PO^O  ^  2*122^2  +  2*233^3  +  ^^®i)  +  ...} 

In  order  to  avoid  singiilarities  for  St  <  0,  the  quadratic  in  (4.5b)  must  be  elliptic  and  of 
opposite  sign  to  ®i.  Since  x^(  is  initiedly  positive,  cf.  (2.1),  it  follows  from  (4.4a)  and  (4.5b) 
that  at  a  first  zero 


*111  >  0  >  *122  >  0  )  *133  >0  >  ®i  <  0  . 


(4.6a,  6,  c,  d) 


The  topology  of  the  characteristics  (4.5),  shown  in  figure  3a,  can  be  compared  to  the 
asymmetric  case  figure  2a,  where  the  separation  characteristic  develops  a  cusp  at  St  =  0. 
In  this  case,  the  separation  characteristic  is  constrained  by  symmetry  to  remain  straight. 

Appropriate  local  scalings  near  separation  can  be  found  using  arguments  similar  to 
those  of  the  previous  section: 


I,  =|6t|i/92ii  =|6t|M2(^"  +  l)>i;  , 
h  =  |^t|Mo/32i2  =  |6t|Mo^2(^’  +  l)’i;  , 

X  =  \St\^^a/3lX  =  \St\^^  001(2^  -h  l)*^*  , 

Y  ¥•  _ 

y=  - r -  =  - 7 - :: - T  »  2=  6t  Z  , 

\st\h^0,  \St\h02{Z^ +1)"  01 

1—  /I -2—  —  \\  P’ 

<1  =  3X111  ,  7  =  (32^3*122*111)  - JJ—  , 

POj-ho* 

(23*122*111)’  \23*111/  V  *111/ 

The  continuity  integral  becomes 

V-  -  t 

Jo  JL*{2X-  -31*  -  L-^)  Jl:  JL-{2X‘  -  3L*  -  L'^) 


where 

iaA'*)  =  /(.V*)  . 

This  can  be  written  as  ein  elliptic  integral  similar  to  (3.16), 


r(£;,.Y-)..  ?f('|m)±  ifMm)  . 


(4.7a) 

(4.76) 

(4.7c) 

(4.7ii,e) 

(‘1.7/,j) 

(4.7h,i,;) 


(4.8a) 

(4.86) 


(4.9a) 
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where 


A(X*)  =  (3(15^  +  +  3))  ^  ,  m(X*)  =  -  - 


1  31;^  +  6 

2  4A2 


^(i;,  J:-)  =  2arctan  +3) 


(4.96,  c) 


(4.9(i) 


Note  that  instead  of  using  LJ  as  the  independent  variable,  there  is  an  advantage  in  using 
XJ,  as  given  implicitly  by  the  relation 


=  {l;^  +  i)i/(xv(x;’  +  i)M  , 


since  at  the  symmetry  line  the  solution  is  regular  in  terms  of  XJ: 

. .  ...  2  /TT 


(4.9e) 


r-(X;.0)~^(^^-arctanX;)  . 


(4.10) 


Contours  of  the  boundary-layer  thickness  F"*'  in  the  X,Z  plane  are  shown  in  figure 
3b.  The  asymptotic  relations  for  large  |.Y|  and  \Z\^  corresponding  to  (3.18)  and  (3.19),  are 

,  (4.n) 


(4.12) 


The  velocity  components  and  density  in  the  neighbourhood  of  the  stationary  point 
are  given  by 

i  ~ -|6<r>|o/?|(Z2  +  1)U;  ,  (4.13a) 


(i,p)  ~  (z,,p,)  -t-  l^t|  ’  (  (2  2,P2)/3o/32(^  +  1)  ’^2  +  (^3,P3) 


■ 


(4.136) 


Hence  XJ  can  be  interpreted  as  the  velocity  component  towards  the  symmetry  plane.  The 
scaled  velocity  profile,  — XJ,  is  illustrated  in  figure  3c  at  a  number  of  X*  stations,  while 
contours  of  XJ,  and  the  corresponding  vorticity  component,  dL\/dY',  are  illustrated  in 
figures  3d  and  3e  respectively.  In  cross-sections  of  constant  z,  the  variations  in  velocity 
parallel  to  the  symmetry  plane  are  proportional  to  XJ.  XJ  velocity-profiles  ^^•e  given  in 
figure  3f,  while  figures  3g  and  3h  illustrate  contours  of  XJ  and  the  vorticity  dL\/dY*. 
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Close  to  the  wall,  i.e.  as  K  — >  0, 

- - ~  “  (^2,^2);^  Y  J  ’  (4.14a, 6) 

which  match  to  a  regiilar  vorticity  layer  of  similsir  form  to  (3.24).  Similarly  a  match  can  be 
achieved  with  a  separating  layer  governed  by  a  Prandtl  transformation  above  the  central 
inviscid  region. 

Figure  3i  shows  the  characteristics,  i.e.  lines  of  constant  x  and  z,  for  6t  >  0.  Inte¬ 
gration  of  the  continuity  equation  yields  a  singularity  over  a  segment  —  1  <  Z  <  1  of  the 
symmetry  line  X  =  0.  Since  neither  the  singularity  nor  any  particles  on  the  symmetry  line 
leave  the  symmetry  line,  the  region  of  inaccessibility  remains  restricted  to  the  symmetry 
line. 

A  special  case  occurs  for  separation  at  the  intersection  of  two  symmetry  lines,  such 
as  at  the  apex  of  an  ellipsoid.  In  that  case,  in  addition  to  the  symmetry  in  a;  is  an  even 
function  of  ^  and  z  an  odd  one,  and  the  transformations  of  the  Lagrangian  coordinate 
system  (4.2)  and  (4.3)  become  trivial.  No  changes  in  the  leading  order  singularity  structure 
occur,  since  it  was  already  symmetric  in  2-direction,  even  though  this  condition  was  not 
imposed.  However,  the  velocity  parallel  to  the  symmetry  plane  must  be  antisymmetric, 
and  the  density  symmetric  (cf.  (4.13b): 

,  p^p,  +  ]8t]^p,/3„02{Z^  +  .  (4.15a,fe) 

Pi  ^3 

In  the  case  of  two-dimensionality,  where  x  is  independent  of  C,  the  coefficient  /?]  van¬ 
ishes  as  in  the  previous  subsection,  suppressing  the  decay  of  the  boundary-layer  thickness 
with  2.  The  flow  on  the  symmetry  line  can  then  be  written  as  a  one-dimensional  problem, 
and  wsis  studied  from  an  Eulerian  standpoint  by  Banks  and  Zaturska  (1979)  and  Simpson 
&  Stewartson  (1982a,b).  In  part  2,  Van  Dommelen  (1989)  uses  this  flow  to  verify  the 
Lagrangian  analysis  numerically  to  high  accuracy.  Favourable  numerical  comparisons  with 
the  singularity  structure  away  from  the  symmetry  line  have  been  obtained  by  Lam  (1988) 
for  starting  flow  through  a  circular  pipe. 

The  existence  of  this  singularity  has  also  been  reported  by  Stern  &  Paldor  (1983), 
Russell  &  Landahl  (1984)  and  Stuart  (1987)  while  studying  inviscid  models  for  the  growth 
of  large  amplitude  disturbances  in  boundary  layers.  In  fact,  because  unsteady  separation  is 
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primarily  inviscid  in  its  final  stages,  an  alternative  approach  to  that  above  would  be  to  solve 
the  inviscid  version  of  (2.3)  exactly,  and  then  to  examine  the  possible  singularities  of  the 
solutions  (see  also  Van  Dommelen  (1981)  for  the  two-dimensional  singularity).  Note  that 
although  Stuart’s  (1987)  exact,  inviscid,  symmetry  line  solutions  do  not  include  a  paredlel 
flow  in  the  ^-direction,  our  results  are  in  agreement  since  the  details  of  the  singularity 
structure  are  independent  of  i(^,  ,t). 

As  in  section  3  the  above  singular  solution  will  not  remain  valid  for  sufficiently  small 
j^t|  because  previously  neglected  pressure  gradients  will  become  important  (cf.  the  inter¬ 
active  problem  for  the  two-dimensional  singulau-ity  formulated  by  Elliott  et  al.  (1983)). 
Further,  because  the  velocity  towards  the  separation  line  is  much  smadler  in  the  upper  and 
lower  vorticity  layers  than  in  the  central  layer,  it  is  in  the  vorticity  layers  that  the  effect  of 
the  pressure  gradient  will  be  felt  first.  However,  it  is  the  central  layer  which  is  responsible 
for  the  growth  in  boundary-layer  thickness;  thus  it  appears  that  the  first  asymptotic  rescal¬ 
ing  does  not  lead  to  an  ‘interactive’  effect  to  smooth  out  the  above  singularity.  Instead,  the 
singularity  continues  to  be  driven  by  the  flow  in  the  central  layer,  while  significant  changes 
occur  in  the  upper  and  lower  layers.  Similar  arguments  seem  to  hold  for  the  singularities 
in  (c)  auid  (d)  below. 

(c)  Axisymmetric  boundary-layer  flow  with  swirl 

In  axi-symmetric  flow,  the  flow  geometry  does  not  depend  on  ^  and  (  individually, 
but  only  on  the  Lagrangian  distance, 

0=\/e-i-  c  ,  (4.16) 

from  the  axis  ^  =  C  =  0.  The  displacement  of  rings  of  particles  g  =  rj  =  t  =  constant  from 
their  original  position  must  remain  restricted  to  a  change  in  physical  distance, 

r  =  \/ ,  (4.17) 

from  the  axis,  a  rotation  around  the  axis,  and  a  shift  in  vertical  position.  Hence  according 
to  the  theory  of  orthogonal  matrices,  the  solution  must  be  of  the  form 

I  =  c(p%rj,t)^  4- 3(e^,Tj,e)(  ,  (4.18a) 

2:  =  -3(p^r7,<)4  +  c(p^TJ,t)(;  ,  (4.186) 
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where,  because  of  the  assumption  of  regular  x  and  z, 


c  ~  -I- |x,^{{{0,T;,0,t)e^  +  ...  ,  (4.19a) 

a  ~  X, ((0,77,0, t)  +  |x.f(f(0,T},0,t)p’ +  ...  .  (4.196) 

In  terms  of  c  and  3  the  physical  distance  from  the  axis  is  given  by 

=  (c^  4-a^)p^  .  (4.20) 

The  Jacobian  J  in  (2.2)  can  be  written  in  terms  of  g  and  r  as 

J  =  i'^^),e^y,v ■  (4.21) 

Thus  separation  occurs  at  a  stationary  point  for  and  from  (4.20)  and  (4.21)  it 

occurs  on  the  axis  when 

c(0,  t;,  ,  t, )  =  3(0,  t;, ,  t , )  =  0  .  (4.22a,  6) 


A  rotation  of  the  Lagrangian  coordinate  system  to  diagonalize  the  second  order  deriva¬ 
tives  of  X  is  not  advauitageous  here,  since  the  axial  symmetry  would  be  lost.  Instead  we 
rotate  the  coordinate  system  around  the  symmetry  axis. 


-  iCuC 

\/xj2  +  X23 


^12^  -I-  g23C 
>/xj2  -I-  Xjj 


(4.23a,  6,c) 


to  eliminate  the  X12 -derivative,  followed  by  the  shearing  transformation 


li=ki  ,  l2=k  +  ~^{^+kl)  +  Stp- 

0X23  X23 


I3  =  ^3 


(4.24a,  6,  c) 


to  eliminate  X333  and  X3. 

The  characteristics  of  the  Jacobian  are  lines  of  constant  distance  r  from  the  axis.  If 
ta  is  the  first  time  that  a  singularity  forms  then  xmij  must  be  negative,  or  for  a  suitable 
choice  of  the  positive  Zj -direction 


xiii  ^  0  ,  X]  <  0 


(4.25a,  6) 


The  characteristic  lines  of  constant  r  in  the  £»,/2 -plane  appear  as  sketched  in  figure  4a. 


Appropriate  local  scalings  are 


5  ,  h  =  ,  T  =  ,  y  =  ,  (4.26a, 6, c,d) 

a  =  j®in  ,  7  =  |®in(®L)*  /’V  >  /^o  =  .  (4.26e, /.^i, /i) 

2(2:53) »  V  ®ni/ 

leading  to  a  continuity  integr^d 


,  r^°  2dP*  I  2dP* 

-(3P*  4-P-3)2  ip.  -  (3P* 


where 


Po*(i2*)  = /(/Z*^ 


(3p.  +p.3)2 


This  can  be  written  as  the  elliptic  integral 


(4.27a) 


(4.276) 


(4.28a) 


where 


A(fi')  =  (3(P;'+l)(P;=+3)>)i  ,  =  1  -  ??»l±ii?2l±i5  ,  (4.284,0 

=2arcUn  +3)'  ,  (4.28d) 

and  P*^  is  related  to  LJ  and  R*  through  the  solution  of  the  cubic  equation 


AR*^  =  +  (P*^  +  3)^P*^ 


(4.28e) 


On  the  axis,  (4.28a)  simplifies  to 


F'(£;,0)  =  ~  +  arclan  Z,;)  , 


(4.29) 


while  for  large  i?*,  the  boundary-layer  thickness  asymptotes  to 


y+  ~ 


4a3  /tt  I  \/3\  I  .  ,2 

forl6t|3«r<l 


3i257 


(4.30) 
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The  velocity  components  in  the  radial  and  azimuthal  directions,  and  the  density,  are 


r^-\st\^lccfi\u;  ,  re^\6tr^la^\{w; -xU;),  X  =  ^  ,  (4.31a, 6,c) 

*1 


k  3  _/j3. 


p  =  t>.  +  -  p)  , 


where  ?  =  sgn(x23),  and 


.  _  (P*»  +  3)P*^  3P*2i; 


u:  = 


,  w;  = 


(4.31d) 


(4.32a, 5) 


are  the  symmetric  and  anti-symmetric  velocity  profiles  shown  in  Figure  4b. 

Both  the  radial  and  the  circumferential  velocity  profiles  depend  non-trivially  on  the 
parameter  x-  However  the  magnitude  of  the  velocity. 


q  =  +w'^  =  |^t|  3  r^a/3^\/l  +  x’P*  , 


(4.33) 


does  not;  contours  of  q  are  illustrated  in  figure  4c.  The  vorticity  components  normal  to 
the  velocity  and  parallel  to  it  aire  proportional  to 


^  U;U;y.+W,-W-y. 

’  ’  ’ 


w*n*,  —rr^w 

^  1^0  ^0^0X1  =  P-(P*2  +  1) 


(4.34o,6) 


(4.35o,fr) 


Contours  of  these  quantities  are  plotted  in  figures  4d  and  4e  respectively. 


A  match  with  the  sandwich  layer  adjacent  to  the  wall  is  eigain  possible,  since  as  F  — >  0 

2  3  1 

r  ~  |7/3^XT'y  ,  r9  --  -\yl3^<iry  ,  p  ~  p,  -  pj - -  (4.36a,6,c) 

3  7  y 

Similarly  a  match  can  be  achieved  with  the  upper  separating  layer. 

Figure  4f  shows  the  characteristics  for  6t  >  0.  The  singular  line  is  the  physically 
expanding  circle  R*  =  1,  but  the  region  of  inaccessibility  is  larger  due  to  particles  with 
Lj  ^  0  which  move  radially  outward  from  the  singular  line  at  a  greater  rate. 

On  the  ajcis  itself,  the  continuity  integral  is  particularly  simple: 


r  P^Hodg 

^  Jo 


(4.37) 
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When  this  integral  is  expanded  to  second  order,  a  logarithmic  correction  to  the  y-position 
is  obtained.  This  and  other  terms  were  initially  overlooked  in  Eulerian  analyses  of  the  flow 
on  the  axis  (Bodonyi  &  Stewaxtson  (1977),  Banks  &  Zaturska  (1981)),  but  in  a  Lagrangian 
approach  the  logarithmic  term  follows  naturally  from  the  hypothesis  that  the  solution  for 
the  motion  parallel  to  the  boundary  is  regular.  (A  similar  logarithmic  term  arises  in  the 
symmetric  case  (b)  above,  cf.  part  2).  In  fact,  from  this  hypothesis  alone,  the  complete 
singularity  structure  presented  by  Stewartson,  Simpson  &  Bodonyi  (1982)  can  be  recovered 
by  means  of  a  simple  integration  of  (4.37). 

(d)  Axisymmetric  boundary-layer  collision  without  swirl 

Finally,  we  consider  the  case  of  axially  symmetric  flow  when  there  is  no  rotation  of 
the  flow  about  the  auds.  In  the  absence  of  such  rotation  (4.18)  simplifies  to 

X  =  ,  2  =  c(Q^,T],t)C  •  (4.38) 


No  transformations  of  the  Lagrangian  coordinate  system  axe  needed  in  this  case.  It  follows 
that  if  a  singularity  first  appears  on  the  axis  c  must  vanish.  The  contours  of  constant  r 
are  then  identical  to  those  for  a  symmetric  collision  (figure  3a),  while  the  Taylor  series 
coefficients  satisfy  conditions  (4.6a,b,d). 

In  a  similar  way  to  before  suitable  loczJ  scalings  are 

,  h  =  ,  r  =  ,  y  =  *  (4.39a, 6, c, d) 

1—  1/1—3—  ^a  ^  f  ^111  ^  f 

a  =  3X111  ,7=3  (33^111X122)^  „  ,  00  =  ,/3=(-z —  I  , 

\XU2  /  \  *111/ 

(4.39e,/,3,h) 

leading  to  a  continuity  integral 


V*  r‘°  P*^dP*  ^  P‘^dP* 

do  R*  y/2R*  -  3P*  -  P*3  /p.  /?•  v^2H*  -  3P*  -  P*^ 


where 


Po(/?*)  =  /(/?*) 


(4.40a) 


(4.406) 


This  solution  can  be  ‘reduced’  to  the  form 


1  /  4  -  4n  TT 


v(r;.«-)=p^  ^ 


2 _ TT 


n(n;  -\m)  -  1^) 
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,2  ^  /<P;  sin(v?)  \ 

2  /nPJ  sin(2v?)  \\  , 


where  n(n;  y?!^)  w  the  incomplete  elliptic  integral  of  the  third  kind,  and 

A{Rn  =  (3(p;’  +  1)(P;^  +  3))  i  ,  m(il-)  =  1  -  ,  n(il*)  =  1  -  , 

(4.416,  c,  d) 

V(i; ,  R- )  =  2  atetan  ^  i  ^ (p;’  +  3)  ^  ,  (4.4l€) 


p-(z;;,R-)  =  (i;' \)\i(r-i(l\^  +  i)S)  . 


2  I  1  \  ?  ' 


(4.41/) 


On  the  axis  (4.41a)  simplifies  to 


while  for  large  R*  the  boundary-layer  thickness  asymptotes  to 

The  velocity  and  density  are  given  to  leading  order  by 

r  ~  -|6<|5|a/?^P*  ,  +|6<|'p2do/3Z:;  • 


(4.42) 


(4.43) 


(4.44a,  6) 


Sample  velocity  profiles  are  illustrated  in  figure  5a,  while  contours  of  P‘  and  the  vorticity 
dP"  jdY*  are  given  in  figures  5b  and  5c  respectively.  Again,  a  match  is  possible  with  the 
wail  layer,  since  as  K  —  0, 


3  5  2  i  i  2 '  1 


3±  i.  k 

»  3  y  3 


(4.45a, 6) 


5.  Discussion 


33 


In  this  paper  we  have  shown  that  the  description  of  attached  flow  paist  a  body  using 
the  classical  boundary-layer  equations  can  break  down  after  a  finite  time  due  the  formation 
of  a  local  singularity.  In  a  Lagrangian  description  the  class  of  singularities  is  characterized 
by  a  singular  continuity  equation,  but  a  regular  momentum  equation.  The  evidence  shows 
that  such  singularities  are  both  mathematically  consistent  and  physicedly  relevant  (e.g.  Veui 
Dommelen  &  Shen  1980a,  Van  Dommelen  1987,  1989,  Lam  1988,  Bouard  &  Coutanceau 
1980).  The  precise  structure  of  the  singularity  depends  on  the  symmetry  of  the  flow, 
and  some  of  the  simpler  structures  have  previously  been  partially  or  totally  described  in 
Eulerian  coordinates  by  other  authors.  The  purpose  of  this  paper  is  to  provide  a  unified 
theory  to  facilitate  the  identification  of  singularities  of  the  Lagrangian  type  when  they  do 
occur.  This  seems  especially  relevant  for  the  difficult  problem  of  the  asynunetric  singularity, 
where  the  singular  behaviour  would  have  to  be  deduced  from  a  three-dimensional  unsteady 
computation. 

These  singularities  occur  when  a  fluid  particle  becomes  compressed  in  one  direction 
parallel  to  the  boundary.  Conservation  of  mass  then  implies  that  the  fluid  above  this  fluid 
particle  is  forced  out  of  the  boundary  layer  in  the  form  of  a  detached  vorticity  layer.  A 
common  feature  of  all  the  singularities  is  that  the  typical  length  scale  in  the  direction 
of  compression  is  0(|5t|»).  However,  the  the  strength  of  the  singularity  increases  with 
the  symmetry  of  the  flow;  the  the  boundary  layer  thickness  varies  from  C)(|5<|~«)  for  the 
asymmetric  symmetry  to  0(|5<|"’)  for  the  axisymmetric  singularity  without  swirl. 

Because  the  singularities  take  the  form  of  a  vertical  ejection  of  fluid  from  the  boundary 
layer,  we  believe  that  they  indicate  the  onset  of  separation  as  hypothesized  by  Sears  & 
Telionis  (1975).  While  the  present  singularity  structures  do  at  least  seem  to  describe  the 
initial  genesis  of  the  separating  shear  layer,  within  an  asymptotically  short  time  interactive 
effects  which  are  neglected  in  the  classical  boundary-layer  formulation  must  be  included 
(e.g.  Elliott  et  al.,  1983,  Henkes  &:  Veldman  1987).  At  that  stage  a  new  asymptotic  scaling 
must  be  substituted  into  the  Navier-Stokes  equations  in  order  to  recover  the  correct  large- 
Reynolds-number  limit.  Knowledge  of  the  precise  asymptotic  structure  of  the  singularities 
is  necessary  to  identify  this  new  scaling,  and  one  of  the  contributions  of  this  work  has  been 
to  identify  the  full  structure  of  a  number  of  symmetric  singularities. 

At  first  sight  the  symmetric  singularities  may  appear  less  likely  to  occur  in  problems  of 


practical  importance.  However,  they  have  previously  arisen  in  inviscid  models  of  ‘transition 
to  turbulence’  in  regions  where  symmetric  counter-rotating  longitudinal  vortices  are  forcing 
the  convergence  of  fluid  particles  (Stern  et  al.  (1983),  Russell  et  al.  (1984),  Stuart  (1987)), 
A  three-dimensionaJ  extension  of  the  work  by  Smith  &  Burggraf  (1985),  may  lead  to 
an  asymptotic  description  of  transition  which  accounts  for  viscosity,  where  the  turbulent 
bursts  are  associated  with  local  regions  of  classical  boundary-layer  separation  (symmetric 
or  otherwise). 

This  work  was  presented  in  part  at  the  lUTAM  “Fluid  Mechanics  in  the  Spirit  of  G.I. 
Taylor”  Conference,  April  1986,  Cambridge.  The  authors  acknowledge  financial  support 
from  the  NASA  “Materials  Processing  in  Space”  program,  the  SERC,  the  AFOSR,  and 
ICOMP,  NASA  Lewis. 
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Figure  1.  Structure  of  the  separating  boundary  layer,  illustrating  the  asymptotic  scalings 
in  the  boundsiry -layer  coordinate  system  (schematic). 

Figure  2.  Structure  of  asymmetric  three-dimensional  separation:  (a)  Lagrangian  topol¬ 
ogy  of  vertical  lines  through  the  boundary  layer  near  the  separation  particle; 
(b)  contours  of  the  scaled  boundary -layer  thickness 

in  scaled,  oblique  coordinates;  (c)  possible  actual  appearance  of  contours  of 
boundary-layer  thickness  (schematic);  (d)  contours  of  the  scaled  velocity  L\  = 
0,  ±1,  ±2 ...  in  scaled  coordinates;  (e)  — LJ  velocity-profiles;  (f)  contours  of  the 
scaled  vorticity  dL\  fdY*  =  0,  ±1,  ±2, . . .;  (g)  topology  of  vertical  lines  through 
the  boimdary  layer  for  times  beyond  the  first  singularity. 

Figure  3.  Structure  of  symmetric  three-dimensional  separation:  (a)  Lagrangian  topol¬ 
ogy  of  physically  vertical  lines;  (b)  contours  of  boundary-layer  thickness  Y'^  = 
3j,3,2 1,. . .  ,1;  (c)  — LJ  velocity-profiles;  (d)  contours  of  =  0,  j,l,  1 
(e)  contours  of  dL\/dY*  =  0,  ±1,±2, ...;  (f)  LJ  velocity-profiles  for  flow  par¬ 
allel  to  the  symmetry  plane;  (g)  contours  of  LJ  =  0,  ±1,±2,...;  (h)  contours 
of  dLlfdY*  =  1,2,3,...;  (i)  Lagrangian  topology  of  physically  vertical  lines 
beyond  the  first  angularity. 

Figure  4.  Structure  of  axi-symmetric  separation  with  swirl:  (a)  Lagrangian  topology  of 
physically  vertical  lines;  (b)  the  velocity  profiles  of  the  two  components  —Uq 
and  Wq]  (c)  contours  of  the  scaled  absolute  velocity  P*  =  0,5,1,!^,...;  (d) 
contours  of  the  scsJed  vorticity  component  normal  to  the  flow  velocity  = 
0,  ±1,±2,...;  (e)  contours  of  the  scaled  vorticity  component  parallel  to  the 
velocity  Qp  =0,-1,  -2, . . .;  (f)  Lagrangian  topology  of  physically  vertical  lines 
beyond  the  first  singularity. 

Figure  5.  Structure  of  axi-symmetric  separation  without  swirl:  (a)  P*  velocity-profiles; 

(b)  contours  of  P*  =  0,  1, 1 1,. . (c)  contours  of  dP"  /dY"  =  0,  ±1,  ±2, . . .. 
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